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ABSTRACT 

We consider the problem of consumption of stars by a supermassive black hole (SBH) at the center 
of an axisymmetric galaxy. Inside the SBH sphere of influence, motion of stars in the mean field is 
regular and can be described analytically in terms of three integrals of motion: the energy E, the 
z-component of angular momentum L z , and the secular Hamiltonian H. There exist two classes of 
orbits, tubes and saucers; saucers occupy the low-angular-momentum parts of phase space and their 
fraction is proportional to the degree of flattening of the nucleus. Perturbations due to gravitational 
encounters lead to diffusion of stars in integral space, which can be described using the Fokker-Planck 
equation. We calculate the diffusion coefficients and solve this equation in the two-dimensional phase 
space {L z , H), for various values of the capture radius and the degree of flattening. Capture rates are 
found to be modestly higher than in the spherical case, up to a factor of a few, and most captures 
take place from saucer orbits. We also carry out a set of collisional A-body simulations to confirm 
the predictions of the Fokker-Planck models. We discuss the implications of our results for rates of 
tidal disruption and capture in the Milky Way and external galaxies. 



1. INTRODUCTION 

The study of collisional relaxation in stellar nuclei 
around massive black holes and the associated rates of 
capture has a long history. The pioneering work of Bah- 
call & Wolf (1976) established a quasi-steady-state solu- 
tion for the stellar distribution, now known as a Bahcall- 
Wolf cusp, which has a density p(r) cx r -7 / 4 inside the 
radius of influence r m , defined roughly as the radius en- 
closing a mass in stars equal to the mass M, of the 
hole. Their solution was obtained from the steady-state, 
one-dimensional Fokker-Planck equation describing two- 
body relaxation and energy exchange between stars in 
the (Newtonian) gravitational field of the massive ob- 
ject, and is characterized by zero (or very small) flux of 
stars with respect to energy into the central hole. 

A more refined treatment requires the concept of a 
"loss cone" , the region of phase space corresponding to 
stars with sufficienly low angular momenta to be cap- 
tured by the bla ck hole: L < L\ c , where the capture 
boundary L\ c ss ^2GM m r\ c is determined either by tidal 
disruption or by direct capture, at some radius r\ c (Frank 
& Rees 1976; Lightman & Shapiro 1977). The lat- 
ter paper also introduced the important distinction be- 
tween empty- and full-loss-cone regimes, with the bound- 
ary between them defined as the energy at which the 
typical change of angular momentum in one radial pe- 
riod, \J (AL 2 ) |T rad , is equal to L\ c . Lightman & Shapiro 
(1977) derived the quasi-steady-state rate of consump- 
tion of stars as a function of energy and showed that the 
distribution function depends logarithmically on angu- 
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lar momentum near the loss cone. These authors, and 
subsequently Bahcall & Wolf (1977), included capture 
from the loss cone via an energy- dependent sink term in 
the one-dimensional Fokker-Planck equation for f(E,t). 
While the addition of such a term greatly increases the 
capture rate, it was found to have little effect on the 
form of the density profile at radii r\ c < r < r m . Cohn 
& Kulsrud (1978) solved two-dimensional Fokker-Planck 
equation in (E, L) space and confirmed the results of 
more approximate one-dimensional studies. 

These early studies were targeted toward massive black 
holes in globular clusters. The theory was subsequently 
applied to determine capture rates in galactic nuclei 
(Syer & Ulmer 1999; Magorrian & Tremaine 1999; Wang 
& Merritt 2004). Studies based on the Fokker-Planck 
equation were also verified by other methods such as 
Monte-Carlo models (Shapiro & Marchant 1978; Fre- 
itag & Benz 2002), gaseous models (Amaro-Seoane et 
al. 2004) and direct TV-body simulation (Baumgardt et 
al. 2004; Preto et al. 2004; Komossa & Merritt 2008; 
Brockamp et al. 2011). 

Relaxation times in galactic nuclei are much longer 
than in globular clusters, and in many cases much longer 
than galaxy lifetimes. Partly as a consequence, galac- 
tic nuclei need not be spherically symmetric, and they 
may contain a substantial population of "centrophilic" 
orbits (saucers, pyramids) that dominate the capture 
rate, even in the absence of gravitational encounters 
(Norman & Silk 1983; Gerhard & Binney 1985; Magor- 
rian & Tremaine 1999; Merritt & Poon 2004; Merritt & 
Vasiliev 2011). In the context of collisional loss-cone re- 
population, by far the majority of studies have assumed 
spherical symmetry. This assumption is not crucial for 
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iV-body integrations (except insofar as it can be difficult 
to construct nonspherical initial conditions), but it is an 
important ingredient of Fokker-Planck studies, because it 
ensures the conservation of angular momentum in the un- 
perturbed motion. To date, all Fokker-Planck treatments 
that allowed for non-sphericity (Goodman 1983; Einsel 
& Spurzem 1999; Fiestas et al. 2006; Fiestas & Spurzem 
2010) have assumed axisymmetry and have further re- 
stricted the allowed form of / by writing / = f(E, L z ), 
with L z the component of angular momentum parallel 
to the symmetry axis. Two integrals of motion are not 
sufficient to specify regular motion in the axisymmct- 
ric geometry however, and numerical integrations in ax- 
isymmetric potentials typically reveal a third integral, -Z3 
(with some orbits chaotic). Of course, the reason for the 
neglect of I3 is the absence of knowledge about its func- 
tional form. For mildly flattened systems, I 3 can be ap- 
proximated by L 2 (Saaf 1968), and this approximation 
has been used as a basis for constructing steady-state 
models (e.g. Lupton & Gunn 1987). 

The neglect of the third integral in the axisymmetric 
problem has several important consequences. Instead of 
individual orbits populated by stars having the same val- 
ues of their integrals of motion, one effectively considers 
ensembles of orbits composed of stars with different I3. 
Moreover, by setting / = f(E, L z ), the phase space den- 
sity is forced to be uniform within this ensemble, which 
may lead to unphysical constraints on the possible evolu- 
tion of the system. The diffusion coefficients must also be 
evaluated as if they did not depend on I 3 , or, more cor- 
rectly, are averaged over all possible values of I3. Finally, 
ignoring the third integral precludes the detailed study of 
regular orbits such as saucers, which might be expected 
to dominate the loss rate (Magorrian & Tremaine 1999). 
However, sufficiently close to the black hole, the unper- 
turbed motion is nearly Keplerian, and standard "plane- 
tary" perturbation theory implies the existence of a third 
integral, which can sometimes even be expressed in terms 
of simple functions (Sridhar & Touma 1999). In this ap- 
proximation, the unperturbed stellar orbits are regular 
(integrable) and respect three independent integrals of 
the motion: E, L z and H, where H is the "secular," i.e. 
averaged, Hamiltonian. 

Given an analytic expression for the third integral in 
the vicinity of the black hole, fully three-dimensional 
Fokker-Planck studies become feasible, describing the 
time evolution of / = f(E, L z , H). For the present study, 
however, we choose to concentrate on evolution in the 
two-dimensional subspace (L z , H) with the gravitational 
potential, and the orbital energy E, fixed. Our justifi- 
cation for ignoring changes in E is the same as in many 
previous studies of the loss-cone problem in galactic nu- 
clei (Syer & Ulmer 1999; Magorrian & Tremaine 1999): 
energy relaxation time scales are typically very long in 
nuclei, too long for steady-state configurations like the 
Bahcall-Wolf cusp to be reached. Instead, the depen- 
dence of / on E is inferred from the observed, radial 
density profile via Eddington's formula. Our goal is to 
generalize the well-known, one-dimensional solutions for 
f(L) in the spherical geometry to the two-dimensional 
case f(L z ,H) in the axisymmetric geometry. 

The paper is organized as follows. In §2 we use the av- 
eraging method to demonstrate the existence of a third 
integral of motion inside the supermassive black hole 



(SBH) sphere of influence and we use it to elucidate 
the behavior of orbits: the tube orbits that are generic 
to the axisymmetric geometry, and the saucer orbits 
that inhabit the low-angular-momcntum parts of phase 
space. The Fokker-Planck equation is written down in 
§3, and a scheme for calculating diffusion coefficients is 
presented in the case of three integrals of motion. Follow- 
ing this general treatment, we then restrict our attention 
to weakly flattened systems, which allows some simplifi- 
cation in the computations. We also concentrate on the 
case p(r) cx r~ 3 / 2 , which is both physically reasonable, 
and which results in analytic expressions for many of 
the diffusion coefficients. §4 discusses the proper bound- 
ary conditions for the Fokker-Planck equation, and §5 
is devoted to the solution of the two-dimensional equa- 
tion and comparison with the one-dimensional (spheri- 
cal) case. It turns out that the flux of stars into the SBH 
is enhanced with respect to the spherical case, but by 
a modest factor: at most a factor of a few. In §6 we 
describe direct A^-body simulations designed to test the 
predictions of the Fokker-Planck models; sections 7 and 8 
briefly discuss the role of chaotic orbits beyond the SBH 
influence sphere, and triaxiality of the stellar potential, 
respectively. Finally, in §9 we estimate capture rates in 
realistic galaxy models, using the Fokker-Planck models 
to access the range of parameters not presently accessible 
to iV-body simulations. 

2. MOTION IN AXISYMMETRIC STAR CLUSTERS 
AROUND BLACK HOLES 

Consider a stellar nucleus in which the density varies 
as a power of radius, with exponent 7, and in which the 
equidensity contours are flattened in the direction of the 
short (z) axis; in other words, an oblate system. Let 
p < 1 be the axis ratio, i.e. the ratio of radii along the 
minor and major axes at which densities are equal. In 
the first approximation, the stellar density and potential 
of a flattened system are described by the spherical part 
modified by an I = 2 Legendre polynomial: 
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where < 7 < 2. The total gravitational potential is 
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where M, is the mass of the supermassive black hole 
(SBH) located at x = 0. 

Throughout this section, we restrict attention to mo- 
tion that satisfies 
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where r g is the gravitational radius of the SBH and r m 
its radius of influence; the latter is conventionally defined 
as the radius of a sphere containing a mass in stars equal 
to 2M». The first inequality permits us to ignore rela- 
tivistic corrections to the equations of motion, and the 
second allows us to treat the force from the stars as a 
small perturbation to the inverse-square force from the 
SBH. The effects of relativity are discussed briefly in Ap- 
pendix A, where more precise conditions for the validity 
of the Newtonian approximation are derived. 

Under these conditions, one expects the motion to be 
nearly Keplerian on time scales comparable with the ra- 
dial period, and we can employ the method of averaging 
(Sridhar & Touma 1999; Sambhus & Sridhar 2000): the 
forces acting on a star are averaged over the unperturbed 
motion, with the orbital elements - the "osculating ele- 
ments" - fixed during the averaging. It is convenient to 
describe the motion using the Delaunay variables which 
are action-angle variables in the unperturbed problem: 
I,L,L Z (actions) and w,w,Sl (angles). Here 



I = icirc = y/GM.a = ^== 



(4) 



is the angular momentum of a circular orbit with given 
semimajor axis a or total energy E, L is the magni- 
tude and L z is the z— component of the angular mo- 
mentum (so that cosz = L z /L gives the inclination of 
orbital plane with respect to the x — y plane); w is 
the radial phase (mean anomaly), w is the argument 
of periastron (uj = corresponds to periapsis in x — y 
plane), and ft is the longitude of ascending node. We fur- 
ther define dimensionless angular momentum variables 
as i = L/I e [0,1] and l z = L z /I e [-£,£}. We will 
also have occasion to refer to their squared values, de- 
noted, following Cohn & Kulsrud (1978), as K = I 2 and 

n z = i\. 

These three pairs of canonically-conjugate variables 
evolve according to Hamilton's equations of motion, with 
the Hamiltonian given by 
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and with expressed in terms of the Delaunay vari- 
ables. We assume that these variables - with the excep- 
tion of the radial phase w - are nearly constant over one 
radial period: 

2 ™ 3/2 

The averaged equations of motion can then be defined as 
the equations of motion corresponding to the averaged 
Hamiltonian 
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where the variables {/, L, L z 
aging of over w. 

After the averaging, H is independent of w and there- 
fore I is conserved, as is the semimajor axis a and the 



energy E. On the other hand, H itself is a (new) inte- 
gral of the motion. Finally, L z is conserved due to axial 
symmetry, from which it follows that the motion is in- 
tegrable. Remarkably, even in the (weakly) triaxial case 
there exist three integrals of motion, L z being replaced 
by another conserved quantity (Merritt & Vasiliev 2011). 

Exact expressions for the averaged Hamiltonian are 
given in Appendix B. A good approximation to the av- 
eraged perturbing potential is 
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which is exact for 7 = {0, 1} (for which Q — 5/2,3/2) 
and approximates the true value to within few percent 
in other cases. This Hamiltonian is very similar to the 
averaged Hamiltonian of the hierarchical restricted three- 
body problem (Kozai 1962; Lidov 1962); a detailed com- 
parison is presented in Appendix C. 

Expressed in terms of a dimensionless time r = 
27rf/TM, the equations of motion read 
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The equations for l z and Q are not needed because i z is 
conserved and because nothing important depends on il. 
The time T M is 1 
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where M*(r) = 47rr 3 p Jr (r)/(3 — 7) is approximately the 
mass in stars within radius r. Tm is the time associ- 
ated with precession of the argument of periastron due 
to the spherically-distributed mass (the "mass-precession 
time"). 

In what follows, it will be convenient to replace H by 
%, a linear combination of H and 1Z Z : 
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To obtain the solution, we solve equation (11) for lo 



1 Note that Tm differs by a factor 2(Q — l)/3 from T prcc defined 
in Merritt & Vasiliev (2011). 
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Fig. 1. — Left: Maximum and minimum possible values of TZ (solid and dashed lines, equation 12b) and precession period T proc (dotted line, 
equation 14) for a series of orbits started with ui = n/2 and TZ = 7?.j n i t , for Ti aep = 0.2 and 1Z Z = 0.05. Orbits with 72.j n it < 7J S cp arc saucers, 
others are tubes. Saucer orbits with initial angular momenta 7^i n it and 7?. z 7?. SC p/7?.i n i t are identical (symmetric about 7£g x = ^/1Z sa p1Z z 
line). 

Right: 1Z(t) plotted for three saucer orbits (equation 13 for lZ Bcp = 0.1, 1Z Z = 0.01). Red dashed, green dashed-dotted and blue dash-double- 
dotted lines are for W = —0.048 (close to the fixed-point orbit), H = —0.02 and H = —0.001 (close to scparatrix). Thin solid lines show 
minimum and maximum possible values of 1Z for given 1Z Z (essentially 1Z Z and 7^ S cp), and dotted line shows 7^fi x = ^/1Z B cpTZz ~ 0.032. 
Another dotted line shows the capture boundary at 1Zi c = 0.015. 



and substitute the result into the first of equations (9): 
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It is clear that 7Z is allowed to oscillate between 7Z m i n 
and lZ max . We thus have two classes of orbit, depending 
on the relation between H2 and IZ3, which are separated 
by the line 7-1 = 0. For % > 0, the orbit is an ordinary 
short-axis tube (SAT), and Hz > IZ2, while in the oppo- 
site case (H < 0) the orbit is a saucer. (A similar class 
of orbits was described by Lees & Schwarzschild (1992) 
for a 7 = 2 scale-free potential; see Appendix D for more 
discussion) . The main distinction is that for saucer orbits 
the angle u) librates around tt/2, which means the apoap- 
sis always lies above (or below) the x — y plane; the orbit 
resembles a conical saucer with inner hole (e.g. Samb- 
hus & Sridhar 2000, Figure 7). For tubes, conversely, ui 
steadily decreases. It is easy to show that saucer orbits 
appear only for TZ Z < TZ S ep (and they only exist in the 
oblate, not prolate case), hence the label "separatrix" . 
The solution to equation (12a) can be expressed ex- 



actly in terms of the elliptic cosine: 
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We call the period of full oscillation in 1Z the "precession 
time." It is given by the complete elliptic integral: 



M 




-Tl n 



For orbits not too close to the separatrix, K ps tt/2. 



(14) 
Note 



that, in a time T prcc , w varies by tt for tube orbits and 
by an amount < tt for saucers. 

It is also useful to write down expressions relating 7Z m \ n 
and lZ max . For tube orbits, 



n z 
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It is clear that if e p <C 1 and 7Z scp < K < 1, these 
two values are quite close to each other, justifying the 
practice of approximating the third integral as L 2 . 
For saucer orbits, the relation is simpler: 



In particular, the condition 7\L m i n = 7vl max 
fixed-point orbit, for which uj = tt/2 always, 
shows the properties of a series of orbits with 
1Z Z < 7?. S ep which start with u = tt/2 and with 
initial angular momenta lZi ni t (left panel), and 
evolution of 1Z for several saucer orbits with 
values of % (right panel). 
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Fig. 2. — 1Z(w) (equation 17) for sets of orbits with the same 
IZz and different H. Top panel: 7£scp = 0.1, 7£ z = 0.15, which has 
only tube orbits; bottom panel: 7^ se p = 0.1, 7£ z = 0.01, which has 
saucer orbits appearing as cycles in the lower part of the plot. The 
fixed-point orbit is marked by the cross. In the upper panel, the 
dashed (red) line is at 1Z = 1Z Z and in the lower panel it marks the 
separatrix. 

Equation (11) can be inverted to express 1Z as a func- 
tion of u) (Figure 2), which is more convenient for the 
averaging procedure in the next section: 
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For tubes only the upper root has physical meaning, 
while for saucers both roots are valid as long as sin 2 u) is 
greater than the following threshold: 
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Finally, we outline the complete phase space in (H,TZ Z ) 
coordinates (Figure 3). The boundary between tubes and 
saucers is % = 0, and the other important boundaries are 
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Fig. 3. — Phase space in (TZ Z ,H) coordinates, for 1Z B , 



0.4. 



The region < H < 1 — 1Z Z is occupied by tube orbits, and the 
bottom left corner by saucer orbits; the lower boundary (blue) cor- 
responds to fixed-point saucer orbits (equation 19). Shown in black 
are lines of constant minimum angular momentum (equation 20) 
for 1Z lc = (0.05, 0.1, 0.15, 0.2, 0.3, 0.4). A dot on the intersection of 
the fixed-point orbit curve and the line of TZi c = 0.1 has coordinates 
given by (21). 

A star is captured if it passes near periapsis while hav- 
ing 1Z < IZic, where 1Zi c = Lf c /I 2 « 2r lc /a is the absorp- 
tion boundary (n is the distance to SBH at which a star 
is either tidally disrupted or captured directly, and a is 
the semimajor axis). In the axisymmetric system, this 
condition corresponds to !Z m i n < TZ\c, although not ev- 
ery star satisfying this condition is immediately lost (see 
§4). The lines of constant 7\L m ; n = 1Z\ C , which mark the 
SBH capture boundary, are straight lines satisfying 
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for tubes 
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(20) 



This boundary in the saucer region touches the fixed- 
point orbit curve (19) at 
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The fraction of saucer orbits in an isotropic nucleus is 
approximately 0.47?. SG p. 

In what follows, we will often assume that lZ scp <C 1 
(or at least not too large), equivalent to assuming that 
the nuclear flattening is modest. 

3. FOKKER-PLANCK EQUATION 

We begin by outlining the general method for deriving 
the covariant Fokker-Planck equation in arbitrary coor- 
dinates (Rosenbluth et al. 1957; Merritt 2013). 
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The local (position-dependent) Fokkcr-Planck equa- 
tion can be expressed in terms of generalized position- 
and velocity-space coordinates (x a , J a ), a — {1,2,3} as 



d[gf(x a ,j a ,t)] 

dt 



+ - 
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i d 2 



2 dJ a dP 



((AJ a )Gf) + (22) 

({Aj a Ajp)gf) 



(Merritt 2013, equation (5.121)). The position and veloc- 
ity coordinates need not be related to each other in any 
particular way (e.g. they not be canonically conjugate): 
for instance, one can adopt the integrals of motion as the 
velocity-space coordinates. In equation (22), coefficients 
in angled brackets represent average and mean square 
changes of the corresponding velocity variables per unit 
time. Equation (22) is valid under the assumptions of (i) 
local encounters that change only the velocity but not 
the position of a star; (ii) weak interactions, which al- 
lows the collision term to be expanded in powers of A J a 
up to second order. 2 Later in this section we will use 
orbit-averaged form of this equation, which additionally 
requires that the significant changes in J a occur on a 
time scale much longer than the orbital period. 

The quantity Q = a 1 ! 2 is the density of states, with a 
the determinant of the velocity-space metric tensor a a p, 
so that the squared distance between two points whose 
coordinates differ by A«7 is given by ds = 
and the number of stars in the phase-space volume 
d 3 Jd 3 x is dN = Qfd 3 Jd 3 x (generalization to an arbi- 
trary, non-trivial metric in coordinate space is straight- 
forward). Under a change of coordinates J a — > J p , the 
coefficients in equation (22) transform as 



r) 1 r) 2 W 



(A>AJ iy ) = (AJ Q AJ /3 ) 



dJ" dJ v 

(Merritt 2013, equation (5.120)) and 
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To compute the diffusion coefficients (AJ a ), 
(AJ a AJ l3 ), one needs to know the distribution function 
describing the field stars. Here we assume that the test 
and field stars are drawn from the same /, but as is 
often done, we replace the exact / by an approximation 
when computing the diffusion coefficients. Namely, we 
neglect changes in the diffusion coefficients caused by 
the evolution, and compute them assuming a field star 
distribution of the form /(£), where £ = —E is the 
binding energy per unit mass. In other words, we neglect 
the anisotropy of the field star distribution. Consistency 
with the (spherically-symmetric part of the) mass model 

2 Bar-Or et al. (2012) argue that retaining only the first two 
terms in the expansion may underestimate the probability of large 
changes in orbital parameters, especially on time scales short com- 
pared with the relaxation time. 
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The assumption of isotropy in the field-star distribution 
is obviously inconsistent with the density model being 
flattened; however, for a small degree of flattening it is 
a reasonable approximation. The flattening may be due 
to streaming motions, to an elongated velocity ellipsoid 
or to both; deviations from isotropy due to nuclear flat- 
tening are of order e p -C 1. 

For an isotropic field-star population, the diffusion co- 
efficients expressed in terms of are well known 
(e.g. Merritt 2013, equations 5.23, 5.55) and may be 
transformed to any other coordinates using (23). We 
first adopt the velocity-space variables J a — {£ ,1Z,1Z Z }, 
then later replace 1Z by % which is a true integral of 
the motion. The expressions for the local diffusion co- 
efficients in these coordinates are given in Appendix E, 
equations (El). 

The orbit-averaged Fokkcr-Planck equation is obtained 
by (i) selecting integrals of motion as the velocity-space 
coordinates; (ii) integrating the local Fokker-Planck 
equation over the phase-space volume filled by an or- 
bit, assuming that f(x a ,J a ) is a constant in this re- 
gion (Jeans's theorem); (iii) using the Leibnitz-Reynolds 
transport theorem to exchange the order of integration 
and differentiation. The result is 



d[G av f(J a ,t)] 
dt 



+ 



d 

i d 2 



2 dJ a 8JP 



((AJ Q ) av g av /) + (26) 
((AJ"AJ^) av g av /) 



(Merritt 2013, equation (5.153)), which has the same 
form as the local equation (22), but now / is understood 
to be a function of the J a and of t only, and the diffusion 
coefficients are averaged according to 



{X) av Q a 



(x) g d 3 x. 



(27) 



orbit 



Setting X = 1 in this expression gives the "density of 
states" C? av , which relates / to the number of stars in a 
velocity space volume element d 3 J: dN = fg av d 3 J. 

In the spherically-symmetric case, orbit averaging re- 
duces to a one-dimensional integration with respect to ra- 
dius, J r r+ Xdr/v r , where r_ and r+ are peri- and apoap- 
sis radii for a given E and L; in other words, averaging 
reduces to weighting X in proportion to the time an orbit 
spends near r. In the axisymmetric case, what is tradi- 
tionally done (Goodman 1983; Einsel & Spurzem 1999) is 
to assume that the orbit fills the configuration-space re- 
gion defined by the condition E-$(R, z)-L 2 z /(2R 2 ) > 
where (R, z) are cylindrical coordinates; in other words, 
existence of a third integral is ignored. Then the average 
turns out to be proportional simply to J J X dRdz, the 
integral being taken over this region. 

In our case, the distinction between saucer and tube 
orbits is critical for the loss-cone problem and we do not 
want to mix orbits having different values of the third 
integral. We must perform the averaging taking into ac- 
count the shapes of the orbits as described in the previous 
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section. This is most easily done by adopting Delaunay 
angular variables (w,u,Ct) as configuration-space coor- 
dinates x a . If the corresponding actions (I,L,L Z ) were 
taken as velocity-space coordinates J a , the Jacobian Q of 
this coordinate system would be unity; but since we are 
using a different set of J", Q is determined by equation 
(24) for the coordinates of choice. 

We split this transformation, and the averaging pro- 
cedure, into two steps. Initially we adopt (£,1Z,1Z Z ) as 
generalized velocities J a and carry out the averaging over 
the radial phase angle (mean anomaly) w, obtaining the 
coefficients for the spherical problem. Since these do not 
depend on u> or fi, averaging over these two angles means 
simply multiplying by (2tt) 2 : 



(X) sph g sph = 4tt 2 / dw (X) det 



f 

Jo 



d{I,L,L z } 



d{£,n,n z } 



(28) 

Here all coefficients, including the Jacobian (explicitly 
written as a determinant), are understood to be func- 
tions of (£,ft, ft z ). In the case of a Kcplerian potential 
$(r) = —GM,/r (i.e. neglecting the contribution of the 
stars), and assuming a power- law density profile for the 
stars with index 7 = 3/2, the field star distribution func- 
tion (in our isotropic approximation) is a constant /o 
(equation 25), and the diffusion coefficients can be ex- 
pressed in terms of elementary functions as 



(Aft) sph = .4 



29 - 66ft 
15 ' 



(29a) 
(29b) 



(ATZ z ) sph - (Aft) sph ^ + ^{9ft(ft - 3ft z ) 

+ 29(1 -K)[(K-K Z ) cos 2 lo-K z ]} (29c) 



(AK 2 ) sph =AK 



58 + 32\/ft-90ft 



15 

ft 2 



(29d) 
(29e) 



<Aft 2 ) sph = (Aft 2 > sp h( ^) +—(K-K Z ) 



2A 



15 



ft 



n 

(A£ATZ) sph = A£ 



x -^[9TZ + 29(1 -ft) cos 2 lj], (29f) 



32\/ft(l - Vll) 



15 

(A£Aft z ) sph = (A£Aft) sph ^, 



(AKAK z ) sph = (A1l 2 } sph 



ft 



where 



A= 167T 2 G 2 TO 2 lnA/ . 



(29g) 
(29h) 
(29i) 

(30) 



Here to* is the mass of a field star (scatterer) and m is 
the mass of a test star (whose evolution is described by 
the Fokker-Planck equation). We set to* = to in what 
follows. 



The density of states G sp h is split into two factors: 



Osph = 



Ge 



2^/WR z 



_ V2*3 {GM.Y 
g £ = £575 ( 31 ) 



so that J Q dlZ J Q dlZ z £/ sp h = Ge is the phase space vol- 
ume element that transforms / to the number density of 
stars M in the spherical geometry: 



Af{£,1l)d£dK = g £ {£)f(£,TZ)d£dTZ. 



(32) 



These coefficients are based on the usual approxima- 
tion of uncorrelated two-body encounters. The effect of 
resonant relaxation (Rauch & Tremaine 1996) would be 
to enhance the diffusion coefficients for ft and 1Z Z ; we 
defer a discussion of this until §6. 

Finally, the diffusion coefficients are expressed in terms 
of {£, %, 1Z Z } by substituting ft(ft, TZ z ,u>) from equation 
(11), transforming the above coefficients under this sub- 
stitution according to equation (23), and averaging over 
the argument of periastron (i.e. the precession phase) u: 



(X) av G av — — 

7T 



<9ft 

duj (X) S ph Gs P h- (33) 



We have explicitly written the Jacobian of transforma- 
tion (24), and the coefficients with tildes are transformed 
from (X) sp h using (23). This is possible because the last 
transformation does not depend on w or f2. The limits 
of integration in (33) are between and it for tubes and 
between u> m i n and ir — w m i n for saucers (equation 18). 

The averaging must be done numerically, however, we 
can derive asymptotic expressions for the diffusion coef- 
ficients in the large- and small-ft regimes. In the former 
case, when ft ^> lZ scp , ft ~ const = H + 1Z Z , so we 
trivially get (AH 2 ) = (Aft 2 ) - 2(AftAft z ) + (Aft 2 ), 
(A^Aft z ) = (Aft 2 ) - (AftAft,). 

In the limit of small ft, i.e. close to the capture bound- 
ary, the asymptotic behavior of the coefficients is differ- 
ent in the tube and saucer regions. It turns out that for 
tubes the most important diffusion coefficient is (AH 2 ), 
while for saucers it is (Aft 2 ) (since the capture bound- 
ary for saucer orbits is almost parallel to the "H axis, the 
flux perpendicular to the boundary is determined mainly 
by the latter coefficient). In the case of tube orbits with 
H + ft z <C ftsep, the asymptotic expressions are 



(Awa>RM 58tt S p 26 Ge 



(34) 



where 



b= 21n2 - -In — 
2 ft, 



sop 



For saucer orbits with 1Z Z <C ft S c pj the asymptotic ex- 
pressions are 



(Aft 2 ) sa 1.5.4 1 + 



H 



ft, 



sop 



ft, 



Ge 



(35) 

Finally, it is convenient to cast the orbit-averaged 
Fokker-Planck equation (26) into flux-conservative form: 



d(G m f) 
dt 



3 W a 3 f 



8J a 
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The drift and diffusion coefficients arc 



V a = -Q w (AJ a ) 



(37a) 
(37b) 



This form is convenient from both from a computational 
and a physical point of view, since (i) it manifestly con- 
serves the total number of stars, and (ii) drift coefficients 
T>tz and T>n z must be zero, which may be demonstrated 
explicitly, or inferred from the natural requirement that 
(in the absence of capture) collisions tend to isotropize 
/, i.e. in the steady state / should not depend on TZ or 

n z . 

4. BOUNDARY CONDITIONS 

As is well known, the steady-state loss rate in the 
spherical geometry can only be inferred from the orbit- 
averaged Fokker-Planck equation if a certain condition 
is satisfied: the change in angular momentum over one 
radial period due to encounters must be small for or- 
bits near the loss cone. Far enough from the SBH, 
this condition is violated (the "pinhole" or "full-loss- 
cone" regime) and the orbit-averaged approximation is 
not valid. This problem is dealt with by returning to 
the local (r, v) Fokker-Planck equation and allowing the 
phase-space density to vary along orbits, from zero at the 
edge of the capture sphere to some finite value at apoap- 
sis (Cohn & Kulsrud 1978). The result is a boundary 
condition for the Fokker-Planck equation that specifies 
the value and slope of j(TZ) at TZ = TZ\ C in terms of its 
average value at the (fixed) energy £ . 

No such exact analysis has ever been carried out in non- 
spherical (axisymmetric, triaxial) geometries, nor will we 
do so here. Instead we will be satisfied with a more 
heuristic derivation of the relation between / and T at 
the loss cone boundary. Our analysis will be similar in 
spirit to that of Magorrian & Tremaine (1999), but, as 
we will see, with some important differences. 

4.1. The spherical "problem 

We begin by reviewing the boundary conditions in the 
spherical geometry. The evolution equation in terms of £ 
and TZ can be derived by integrating equation (36) over 
1Z Z from to 1Z. We introduce two, two-dimensional 
quantities: the number density of stars, 

Af(£, TZ) d£dlZ = / dlZ z G sph (£,lZ,TZ z ) f{£,1Z,1Z z )d£dTZ 
Jo 

= 4n 2 T lad (£)L 2 circ (£) f(£,TZ) d£dTZ 
= g £ (£)f(£,lZ)d£dlZ, (38) 
and the flux per unit energy in the 7?.-direction: 

Fk(£, n)d£=- j n dTZ z g sph d£ 

(ATZ 2 ) dN{£,TZ) 



dTZ 



d£. 



(39) 



The relation between Tji and the fluxes appearing in 
equation (36) is that the former is the integral, over the 
entire loss cone boundary, of the component of F a nor- 
mal to that boundary, in the subspace £ = const. In the 



spherical case, the capture boundary in the 1Z — 1Z Z plane 
is defined by 1Z = TZ\ C , and so T-r. = J Q K dTZ z F n . 

As emphasized by Frank & Rees (1976), in the spher- 
ical geometry, changes in angular momentum are ex- 
pected to dominate the capture rate. Setting £ = con- 
stant, we consider one-dimensional diffusion in 1Z, which 
obeys 



dAf(K,t) _ dT; 



n 



dt dTZ 
1 d ( o dM\ 

~ 2 87z{ {An) mz)- 



(40) 



The final expression uses the result that for small 1Z, 
(ATZ) = \(d/dTZ)(ATZ 2 ). Also in this limit, (Aft 2 ) cx 1Z, 
and we can write 



V{£) = lim 



(ATZ 2 ) 
2TZ 



(41) 



an inverse, orbit-averaged relaxation time. 

If we adopt the approximation (ATZ 2 )/ 2 — T>TZ over 
the entire range of TZ, then equation (40) is mathemath- 
ically equivalent to the heat conduction equation in a 
circular domain (Ozisik 1993). The steady-state solu- 
tion is N cx \nR + const. The natural boundary condi- 
tion at TZ = 1 is T-r. = 0, which in this approximation 
translates to dN /dTZ = (in the more exact treatment, 
(ATZ 2 ) itself tends to zero at TZ = 1). The boundary at 
TZ = TZ\ C is responsible for capture, however, simply set- 
ting N{TZ\ C ) = is not always appropriate. As discussed 
by Lightman & Shapiro (1977), there are two regimes 
characterizing the behavior of Af near loss cone bound- 
ary, depending on the ratio q between the radial period 
and the time to random- walk out of the loss cone (that 
is, to change TZ by order of TZ\ C ): 



rp _ K = TZ ± 

™ ~ (ATZ 2 )/2 V ' 



(42) 



The case q <C 1 is the "empty loss cone" regime, since a 
star scattered into the region TZ < TZ\ C is swallowed much 
faster than encounters can scatter it back out. The op- 
posite situation, q ^> 1, is the "full- loss-cone" regime, 
because diffusion is so rapid that even capture of stars 
(near periapsis) does not substantially diminish the pop- 
ulation of loss cone orbits away from periapsis. 

It turns out that much of the flux into the SBH comes 
from stars at energies where q w 1, so the behavior of the 
solution at the transition between the two regimes is im- 
portant. Let Mc = N{TZ\c), and define F\ c = — Jrc(7£i c ), 
the flux through the boundary (i.e. the number of stars 
captured per unit time per unit energy). In a steady 
state, T-ji is independent of TZ near TZ\ C . We can express 
the inner boundary condition in a general way as 



•Fie 



a 



(43) 



where all quantities arc understood to be functions of 
£, or equivalcntly of q(£). After usinq equations (38) 
and (39) to express M and T-n in terms of / and 
df/dTZ, equation (43) is seen to be a boundary condi- 
tion of the Robin type (linear combination of function 
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and its derivative) (Eriksson et al. 1996). The dimen- 
sionless coefficient a can be derived by returning to the 
local (non-orbit-averaged) Fokker-Planck equation and 
determining how / varies with radial phase assuming 
/ = at periapsis (Baldwin, Cordey & Watson 1972). 
Cohn & Kulsrud (1978) produced a numerical solution 
in the spherical geometry and proposed the approxima- 
tion a — 0.186g + 0.824y / <7 for q < 1 and a = q for q > 1. 
An exact solution exists to the same set of equations 
solved by Cohn & Kulsrud: 



~<Kn.ll A 



(44) 



(Mcrritt 2013), where a m are consecutive zeros of the 
Bessel function J . Equation (44) is unwieldy; a good 
approximation is 



«W = (9 4 + <z 2 ) 1/4 

which has the asymptotic forms 

'^/q + q 5 / 2 /4 if<7«l 
q+l/(4q) if«»l. 



(45) 



(46) 



These expressions differ most strongly near q = l, where 
the exact solution (44) gives 1.195. Cohn & Kulsrud's ap- 
proximation is a(l) = 1 while equation (45) gives 1.189. 

In terms of a, the variation of AT with 1Z near the loss 
cone boundary is 



lZ>TZo = 1Z\ C exp(— a). 



(47a) 
(47b) 

Here IZo < 1Z\ C plays the role of the effective absorbing 
boundary at which 7f = 0. 

It is convenient to introduce the "draining rate" of a 
uniformly-populated loss cone: 



*^*drain — 



77\ C 1Z\ C Sf/lc^lc 



T, 



rad 



T, 



rad 



= A^LUxc (48) 



This is the capture rate that results if the following two 
conditions are satisfied: (i) the distribution function is 
constant inside the loss cone, with value /i c ; (ii) the vol- 
ume of the loss cone (per unit energy), Gs1Z\ c , is emptied 
every radial period. As is apparent from its definition, 
•^drain does not depend on the diffusion coefficient T>. 
We can rewrite the boundary condition (43) in terms 

Of J~ drain 

(49) 



T\c = — J^dr 

a 



In the full-loss-cone regime, q w a and the capture rate 
is Fie ss .Fdrain; otherwise q/a < 1, reflecting the fact 
that the phase space density decreases to zero at some 
K = 1Z <Ki c (47b). 

Dchnc 77(£) to be the integral of 7f(£,1Z) over angular 
momenta: 



77(£) = f 77{£,iz)diz. 



(50) 



Roughly speaking, AT is the quantity that would be in- 
ferred from an observed radial density profile (via Ed- 
dington's formula, say). If we extrapolate the logarith- 
mic solution (47) to all 1Z, we can relate N(1Z) and the 



capture rate T\ c to 77: 

N(1Z)k — -A, (51a) 
a + ln(l/7Ci c ) - 1 



v 77 



a + ln(l/n h 



1 



(51b) 



In the full- loss-cone regime, q 3> 1, the distribution func- 
tion is almost isotropic (A/ (1Z) w 77) , and 



- rad 



(52) 



That is: the full volume of the loss cone is consumed 
every radial period, and ^drain and F\ c are equivalent. 
This also means that the exact value of the diffusion co- 
efficient or even the very process responsible for shuffling 
orbits in 1Z does not affect the capture rate, as long as 
it is efficient enough to keep the loss cone full. In the 
opposite limit of q <C 1, 



•Fie -> 



77TZ h 



ln(l/ftl c ) Trad 



(53) 



Now the capture rate is limited by diffusion from larger 
1Z, that is, by the gradient of 77 near 7Zi c . The distribu- 
tion function is depleted at small 1Z. 

Here we note a distinction that will be important when 
discussing the axisymmetric problem. Equation (48) ex- 
presses .Fdrain in terms of the value of 7f at the loss cone 
boundary. In the spherical case, the full-loss-cone bound- 
ary condition (a « 5 > 1) implies 77\ c w Af (equation 
51a). This is because the same process - gravitational 
scattering - is responsible both for populating loss-cone 
orbits uniformly with respect to radial phase and for driv- 
ing the global shape of f(lZ) towards isotropy. As we will 
see, the same is not necessarily true in the axisymmetric 
geometry, because these two actions are driven by differ- 
ent processes: the latter is always attributed to two-body 
relaxation (scattering), but the former may also be driven 
by regular precession. In what follows, we use the terms 
"empty" and "full" loss-cone regimes to distinguish be- 



tween the cases when flux F\ c <C ^d 



and w Ta 



respectively, whatever the global shape of the solution. 

In the spherical problem, the transition between the 
two regimes is naturally defined as a rts q = lnl/7\Li c 
(so that expressions (52) and (53) are equal). Some- 
times another definition is used, based on the require- 
ment that the draining rate equals the repopulation rate 
from nearby regions in phase space, implying q = 1. In 
§6 we denote the energy of the former transition as £ g i bai 
and the latter as £i OC ai> m reference to the fact that these 
conditions are based either on the global shape of the so- 
lution or on its local properties near the capture bound- 
ary. 

Returning to the time-dependent equation (40), if 
(A1Z 2 ) = 2VTZ, an analytic solution exists in terms 
of Bessel functions (Milosavljevic & Merritt 2003; Mer- 
ritt & Wang 2005). If we take 77(11) = 0(1Z - 1Z Xc ) 
as the initial condition, where G is the Heaviside step 
function, then the logarithmic profile is established after 
At w 10 _2 I? _1 , and the flux F\ c , after the initial tran- 
sient, is well described by the quasi-steady-state value 
(51b). Numerical solutions to equation (40) without the 
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simplifying assumption (A7Z 2 )/7Z — const are found to 
match the analytical solution very well (to within a few 
percent). In § 6 we will refer to both the quasi-stationary 
flux value and the time-dependent solution. 

One should keep in mind that the capture rate in 
the time-dependent case may be substantially higher 
than the stationary flux for t <C 1/2? , at least in the 
cmpty-loss-cone regime; sometimes the ratio can be more 
than an order of magnitude (e.g. Milosavljevic & Mer- 
ritt 2003). This, however, depends critically on the de- 
tails of the initial distribution. If instead of a step- 
function at TZ\ C one starts from a profile with a larger 
area 7^-dcpi 3> TZ\ C where the distribution function has 
been depleted, for example, as a result of a binary SMBH 
scattering away stars with periapsides smaller than the 
binary separation, then initially the capture rate is, con- 
versely, much lower than the stationary flux (Merritt & 
Wang 2005). 

4.2. The axisymmetric problem 

We first summarize the various time scales in the ax- 
isymmetric geometry. The three times defined above that 
characterize motion in the smooth potential satisfy 



Trad <C Tm < T p 



(54) 



where T ra d (equation 6) is the radial (Kcplerian) period, 
Tm (equation 10) is the approximate time for uj to change 
by 2ir due to the spherically-distributed mass, and T prcc 
(equation 14) is the oscillation time for 1Z due to torquing 
from the axisymmetric component of the potential. The 
latter inequality is strongest for saucer orbits that are 
near the separatrix and which precess very slowly (Figure 

1). 

The two-body relaxation time can be estimated from 
the diffusion coefficients by taking the inverse of the com- 
mon dimensional factor A, equation (30): 



T rc i = A 1 = Tt 



m 



M. 4\/2 
to* 97rlnA 



(55) 



This time is smaller by a factor 0.87 than the more stan- 
dard definition of the relaxation time in terms of local 
density and velocity dispersion (cf. Merritt (2013), equa- 
tion 5.61). Clearly T re] > T M . 

Magorrian & Tremaine (1999) used the term "loss 
wedge" to describe the set of orbits which can be cap- 
tured by the SBH in the absense of relaxation, i.e. if their 
angular momentum falls below 1Z\ C at some phase of pre- 
cession cycle, or, equivalently, if TZ m { n < 1Z\ C . The name 
reflects the fact that this region is elongated in the H di- 
rection much more than in 1Z Z (Figure 3): its boundary is 
defined by setting TZ m i n (H,lZ z ) = 1Z\ C in equation (20). 
In what follows, we define /i w as the value of the dis- 
tribution function at the loss wedge boundary, which we 
approximate to be constant throughout the loss wedge, 
and .Fi w as the capture rate of stars per unit energy: 



F n >dU 



"Hlc.FPO 



(Ft 



H 

tube 



-^saucer) d1Z z 



(56) 

Here F n and F n * are fluxes defined in equation (36), 
"Hic,fpo is the lowest possible value of % (21), and the 
two terms in the last integral give the contributions to the 



capture rate from the "downward" flux in the % direction 
in the tube region of phase plane and the "upward" flux 
in the saucer region (see Figure 3). Signs are chosen so 
that F\ w is the positive rate of capture. As discussed 
near the end of the previous section, in the tube region 
F u is the dominant flux, while in the saucer region F n * 
dominates. In the case 7Z\ C <C 7Z sep , most of the loss 
wedge lies in the saucer orbit region, and in this region 
the largest contribution to F\ w comes from the flux in 
the 1Z Z direction. 

The number of stars (per unit energy) inside the loss 
wedge is given by 3 

W<i w = J j Q m fdUd1Z z M G £ f lw ^K scv K Xc , (57) 

where the boundaries of the region of integration are 
given by 1Z Z < 1Z\ C and equation (20), and we have taken 
/ to be constant (/i w ) within this region and used the 
asymptotic expression (35) for Q av in the saucer region 
(which gives the main contribution to the integral in the 
case 1Z\ C <C 2£ sep ). The number of stars which are in- 
stantaneously inside the loss cone {TZ < 7Z\ C ) is essen- 
tially the same as in the spherical case: GsfiwTZic- This 
is an expected result: in the axisymmetric case, the ef- 
fective volume of the loss region increases by a factor 
~ V^sep/^ic; but the probability for any star inside 
this loss region having an angular momentum less than 
the capture threshold decreases by the same factor. 

We are now in a position to derive the boundary con- 
dition, i.e., to relate the capture rate F\ w to the value of 
/ at the boundary of the loss wedge, f\ w . An orbit inside 
the loss wedge is captured if its instantaneous value of 
TZ is less than 1Zi c , i.e., if it is in the loss cone, during 
periapsis passage. Here, as in the spherical case, there 
are two possible regimes. If the radial period is short 
compared with the time required for an orbit to pass the 
minimum of its precession cycle while having 1Z < 1Z\ C , 
then every orbit in the loss wedge will be captured in a 
time no longer than one precession period. We call this 
the "empty loss wedge" regime. The rate of consump- 
tion of stars per unit energy, F\ W) is then given by the 
number of stars inside the loss wedge, 7V<i w , divided by 
their lifetime on these orbits, T prcc : 



A/', 



rain.lw 



<lw 



prcc 



(58) 



In the opposite limit, a star that achieves 1Z < 1Z\ C 
while being far from the SBH may precess out of the loss 
cone before reaching periapsis, similar to what happens 
in the full-loss-cone case of the spherical problem. Then 
the capture rate is less than given by the above equation, 
because not all stars in the loss wedge are captured after 
one precession period. It is easy to see that in this case, 
which can be called "full loss wedge regime" , the rate of 
consumption is equivalent to the draining rate of the full 
loss cone (48). In other words, the precession is fast and 
shuffles stars in angular momentum quickly enough that 
the loss cone stays full, hence the capture rate is just 
the instantaneous number of stars inside the loss cone 

3 Note that here JV<i w denotes the integral of the distribution 
function over the loss region, not its value at the boundary as in 
(38). 
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divided by their radial period. By ignoring the effects of 
a finite precessional time, Magorrian & Tremaine (1999) 
were essentially in this regime. 

By analogy with the spherical case, we introduce the 
quantity g ax i separating the two regimes: 4 



rain,lw 

#axi = 



' rad 



drain 



T 



K 



scp 



(59) 



It is easy to see that <7 ax ; 3> 1 at the radius of influ- 
ence, where T rac j ss T prcc : unlike the spherical case, the 
transition from empty- to full-loss-wedge regimes always 
occurs well within the radius of influence, and therefore 
the main contribution to the total capture rate comes 
from the full-loss-wedge regime. Moreover, for most re- 
alistic cases g ax i 3> q for the entire range of radii. Indeed, 
combining expresssions (10, 42, 55, 59) and substituting 
1Z\ C = 2r\ c /a, where a is the orbit semimajor axis, we 
obtain 



q&xi 
q 



\/Tt scv 1Z\ c ^ lZ scp M. 
DTprec 5 In A TO* 



(60) 



This ratio decreases with radius; evaluating it at a = r n 



and substituting n c > 8r g = 8GM./c 2 



GM./a 2 



where a is the velocity dispersion of stars outside r m , 
we may estimate each term in the above expression as 
<7axi/<7 ^ 10~ 2 10 6 10~ 3 > 1. In other words: changes 
in angular momentum near the loss cone boundary are 
determined by the regular precession (q ax i), and not by 
relaxation (q). 

To summarize, the boundary condition in the axisym- 
mctric problem is 



(61) 



and the relation between the flux and the value of / at the 
boundary, expressed in the same way as in the spherical 
problem (43), reads 



Qs /lw = a a xi£> 1 Jiw, 

'q/q&xi if <? ax i < l, 

a if <? axi > 1. 



C^axi — 



(62a) 
(62b) 



The derivation above only gives this relation in "inte- 
grated" form, that is, one coefficient a ax ; for the entire 
H — 1Z Z plane at a given energy. While this is certainly 
an oversimplification, we argue below that it does not 
greatly affect the capture rate. 

The distinction between empty and full loss cones in 
the spherical problem depends on whether Ti c <C .Fdrain 
(equation 49) or not, or equivalently whether a q. In 
the spherical case, a w gmax(l, q^ 1 ^ 2 ), and the transi- 
tion occurs at q s=s 1. In the axisymmetric problem, the 
distinction between empty and full loss wedges is whether 
T\ c <C .Fdrain or not, and the transition is at g ax j = 1. In 
most realistic cases, q ax i » q, although it is not necessar- 
ily true that g ax i > q 1 / 2 . Therefore, the coefficient a ax i 
in the boundary condition (62a) can be both greater or 
less than its spherical counterpart a (equation 45) for 
the same £ and 1Z\ C . In the case q > 1 (full loss cone 
of the spherical problem) there is essentially no differ- 
ence in boundary conditions since g ax i is also > 1. In the 

4 Merritt & Vasiliev (2011) defined an analogous quantity for 
pyramid orbits in the triaxial geometry, their equation (55). 



opposite case (q < 1), a ax i < 1 regardless of the value 
of qaxi, and the capture rate turns out to depend only 
weakly on it, as argued in the next section. In this latter 
case <7 ax i may be both greater or less than unity, i.e. the 
loss wedge may be either empty or full. However, as we 
noted above, the relation between the capture rate and 
/ in the axisymmetric problem depends not only on the 
boundary condition, but also on the structure of the en- 
tire solution, which will be addressed in the next section. 

5. SOLUTION OF THE TWO-DIMENSIONAL 
FOKKER-PLANCK EQUATION 

We are primarily interested here in the capture rate, 
i.e. the flux of stars into the SBH, and not in the evolu- 
tion of the mass distribution (density profile, flattening) 
which we assume to be fixed. The flux is determined 
mainly by diffusion in angular momentum. Accordingly, 
we consider the two-dimensional Fokker-Planck equation 
describing evolution in (H,1Z Z ) and neglect diffusion in 
energy. We note that most of our results are quoted for 
a density profile p oc r~ 3 / 2 which is reasonably close to 
the Bahcall-Wolf stationary solution, p oc r~ 7 / 4 , further 
justifying the neglect of energy evolution. 



5.1. Pr 



studie 



To date, all studies of axisymmetric systems assumed 
that the distribution function depends only on the two 
classical integrals of motion, £ and 1Z Z . Then it is easy 
to show that f(1Z z ) oc \fW z + const for R z « 1. In- 
deed, from (29f) we see that {A1Z 2 } oc 1Z Z for small 1Z Z , 

— 1/2 

and from (31) that the density of states Gav oc 1Z Z ' . 
Then the total capture rate per unit energy is J 7 oc 
Gav (A.1Z 2 .) df/dTZz, and should be independent of 1Z Z , 
which leads to the square- root profile of f(lZ z ). 

Magorrian & Tremaine (1999) used this argument to 
derive the relation between the capture rate F and the 
average value of the distribution function at a given en- 
ergy /, as follows. Start by writing the relation be- 
tween T and the value of / at the loss wedge boundary 
/iw = f(1Z z = TZ\c), which corresponds essentially to the 
full loss wedge regime (equation 62a with g axi = 1). Then 
express the integrated flux in the 7l z direction as 



(63) 



independent of 1Z Z . The numerical factor B is related to 
the "area of the loss wedge" , B w ^TZ m /n, where lZ m is 
the peak angular momentum of saucer orbits and may be 
associated with our definition of 7?. sep . The distribution 
function is then 



OT , / 2 

m z ) = f lw + ^- £ Vn z = flw (i + -V^ 



The average distribution function is 

J(£) = f dn z w(n z )f(n z ) , 

Jo 



(64) 



(65) 



where w(lZ z )dlZ z is the fraction of the phase space vol- 
ume at a given 1Z Z . Magorrian & Tremaine (1999) took 
w = 1/(2-^/7^7); a more correct value is w = — 1. 
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Fig. 4. — (1Z Z ,H) phase plane showing stream lines in a quasi- 
stationary solution to the two-dimensional diffusion problem, for 



0.25 and Tt\ c = 0.003. The capture boundary 7^. m i n ('H, Ti z ) 



(red) stretches from H ~ — T£ scp to H = 1Z\ C almost parallel to the 
ordinate (for an exagerrated close-up refer to Figure 3 where lines 
of constant Ti m i n are shown in solid black). More than one- half of 
the flux lines end up in the saucer region (at H < 0). Fill color 
shows the value of /. 

Using their value, one finds 

the correct expression would contain 1 /3 instead of 1 /2 in 
the brackets. The relation between the steady-state cap- 
ture rate per unit energy and the average (isotropized) 
value of the distribution function is 



Vjf 



MT 



(66) 



where J7 =Q £ f. 

Comparing of this expression with the analogous one 
in the spherical case (51b), we see that when g> 1, the 
capture rate is essentially the same as in the spherical 
case (full- loss-cone regime) , while in the opposite limit it 
is determined by the diffusion coefficient T> and the value 
of lZ scp , rather than by the size of the loss cone TZ\ C . This 
is a consequence of the geometry of loss wedge boundary, 
which stretches in one direction to a fixed fraction lZ scp 
of the phase space. 

5.2. The present study 

We set up an initially uniform (isotropic) distribution 
(/ = const) in the H — 7Z Z plane outside the capture 
boundary, defined by equating 1Z\ C and 7\L m i n found from 
equation(12b) (IZ2 for saucers and IZ3 for tubes); a series 
of isolines of constant 7l m i n is plotted in Figure 3. We 
studied a range of values for both 7\L SGp and 1Z\ C <C lZ sep , 
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Fig. 5. — Steady-state numerical solution of 2d problem for 



0.1, Ti- 



le 



3 X 10~ 



0, together with asymptotical 
profiles at small 1Z Z and large v = H + Ti z . Red short-dashed 
line is f(u) in the tube region, averaged over lines of constant v 
(which roughly correspond to average 7Z for v > TJ. S ep); purple 
long-dashed line is /(Ti z ) in the saucer region, averaged over lines 
of constant 1Z Z ■ Blue dotted line is the approximation (68) for the 
saucer region, and solid cyan line is the approximation (69) for the 
main (tube) region of phase space, which intersects the abscissa 
axis at the effective capture boundary Ti c ff ~ 0.015 (72). These 
two asymptotical profiles do not intersect at (TZ B ep , /sep ) as they 
would if we used the crude approximation for 1Z c ff described in the 
text; the actual value for lZ e ff was computed from the numerical 
2d solution and is roughly twice as higher than the simple approx- 
imation. Green dot-dashed line shows the value of /(Ti z ) for the 
entire phase space, averaged over the second coordinate; it is not 
at all close to square-root profile used by Magorrian & Tremaine 
(1999) and does not tend to zero at 1Z Z = Tii c , since most part of 
phase space at fixed TZ Z is occupied by tube orbits which do not 
come close to the capture boundary. 

as well as various parameters a ax i in the boundary con- 
dition (62b). 

The numerical solution of equation (36) was obtained 
on a non-uniform rectangular grid using two different 
sets of coordinates, defined such that the capture bound- 
aries are parallel to the coordinate axes (see Appendix 
F for details); grid sizes were typically 100 — 300 in each 
dimension and the loss region was resolved by 10-20% 
of the grid cells. We advanced the solutions until time 
T = 1/2? to achieve a steady-state profile, from which we 
could extract the relation between the capture rate and 
the average value of /. 

In the spherical case, the solution is controlled by two 
parameters (aside from T> which scales the time): the 
capture boundary 1Z\ C and the boundary coefficient a (or 
q). As regards the steady-state profile, these two parame- 
ters are not independent, since one can always transform 
to a problem with a' = and TZ[ C = IZo (equation 47b), 
so the family of solution is effectively one-parametric. 
To compare the time-dependent solution of the axisym- 
metric problem to the spherical case, we introduce the 
concept of "equivalent spherical problem", that is, the 
one-dimensional problem with the same coefficient a in 
the boundary condition as a ax i in equation (62b), and 
with some effective capture boundary lZ e g chosen such 
that the time-dependend capture rate closely follows that 
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B « OA^/TZscp (which is ~ 30% larger than the value 
used by Magorrian & Tremaine (1999)). The solution in 
the saucer region is obtained by solving the differential 
equation (67): 



/saucer (^2 ) — /lw~l"" 



IF 

Bvg £ 



n z i 




(68) 



Fig. 6. — Stationary fluxes and effective capture boundaries for a series of 2d problems with a = (empty-loss-cone limit). Top curves 
(boxes) have 7^. se p = 0.25; middle curves (circles) are for 7^sep = 0.1; bottom curves (triangles) have 7^. S cp = 0.03. 

The left panel shows the ratio of the stationary capture rate to the average value of distribution function, T I f. Red dashed curves are 
the capture rate from tube orbits, green dot-dashed curves from saucer orbits, and blue dotted curves are the total capture rate. The blue 
solid curves show the approximation of equation (73). 

The right panel plots the effective capture boundary 7£ e ff normalized to TZ Bcp - Blue dashed lines are results fitted from the 2d solution 
(best-fit logarithmic profile for Ti > 1Z Bep; like the green dot-dashed line in Figure 5); blue solid curves are the approximation of equation 
(72). 

of the axisymmetric problem. Our goal is then to find 
1Z e g as a function of the loss cone size TZ\ C and degree of 
flattening, the latter parametrized by lZ scp - 

In the remainder of this section we present simple an- 
alytical arguments that give a qualitatively correct de- 
scription of the two-dimensional numerical solution of 
the axisymmetric problem, and provide a fitting formula 
for TZ e s. 

Figure 4 shows stream lines of flux and isocontours of 
/ in a quasi-stationary 2d solution for a rather exagger- 
ated value of IZic — 0.003 ~ 10 _2 7?. SC p. Even in this case, 
most of the stream lines end inside the saucer region, 
and that is definitely so for more realistic (smaller) val- 
ues of 1Z\ C . It is also clear that in the saucer region, / 
depends mainly on 1Z Z and is almost independent of the 
second coordinate, which justifies the square-root profile 
of /(72-z) as in equation (64), but only in this region. In 
the tube region, for the greater part of the phase space 
(v = H + 1Z Z > H S ep), the solution is close to that of 
the spherical problem, that is, f(TZ) oc In 1Z + const, with 
1Z m v experiencing only small oscillations. We can build 
an approximate solution by joining the two asymptotic 
forms at 7?. scp . A better description for the saucer region 
accounts for the fact that the flux in the TZ Z direction 
gradually decreases from F at the capture boundary to 
zero at 1Z=, 



which is a somewhat improved form of equation (64). 
The solution in the tube region outside TZ scp is approxi- 
mated by 



/tube (7^) /sep 



(69) 



the coefficient of the logarithmic term gives the same 
flux in the 1Z direction as in the spherical problem, and 
/sc P = /saucer (^sep)- Figure 5 shows the two asymptotic 
expressions along with the actual numerical solution. 

We compute the isotropized value / taking into ac- 
count only the contribution from /t u bc 7 which introduces 
a fractional error of at most TZ S cu- 



f 



/tube(ft) dfc « /, 



scp 



F 



111 - 



1 



- 1 



^scp ■ 



,0 i 

Fn, = / dn-(ATZ 2 z )G a 

" 'Hmin 



df 
dK z 



Fx 11- 



(70) 

Putting all this together and expressing the relation 
between F and N in terms of the coefficient a ax ; (62a) , 
we obtain 



(67) 

T-Lmm{TZz) is the minimum value of H, corresponding to 
the fixed-point saucer orbit (19). The second, approx- 
imate equality in equation (67) is an empirical fit to 
the numerical 2d solution. Using the asymptotic expres- 
sions (35), it is easy to show that for TZ Z <C 7£ sep the 
flux has the form (63), with the numerical coefficient 



F : 



Vjf 



ln(l/ft sop ) - 1 + 2^/K^ P ~/B 



(71) 



Equation (71) can be compared with equation (66) of 
Magorrian & Tremaine (1999): both share the property 
of being independent of 1Z\ C1 replacing it with some ef- 
fective capture boundary 7?. e fT for the empty-loss-cone 
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regime, although this effective value is different. By com- 
paring (71) with the spherical analog (51b), we see that in 
our approximation, TZ e s — lZ sop exp(— l/B) w 0.087^. SC p. 
Figure 6 shows that equation (71) predicts well the flux 
in the numerical steady-state solutions; a better approxi- 
mation to the effective capture boundary and the capture 
rate is 



K eS = Tl sep 0.1 + 0.9 

\ \/ 'v 



scp 



T = 



VN 



a axi + ln(l/7£ eff ) - 1 ' 



(72) 
(73) 



Overall, the capture rates in the axisymmetric geome- 
try are higher than in the spherical case with the same 
boundary condition a = a ax ;, but not by a large factor: 
in the full-loss-cone regime (a 3> 1) they are essentially 
the same, while in the empty-loss-cone regime the effec- 
tive boundary TZ e s is higher than 1Zi c , but since the flux 
depends on it only logarithmically, the difference is not 
likely to be more than a factor of a few. 

Of course, more relevant is a comparison that takes 
into account that a in the spherical case may be different 
from a ax i for the same values of £ and TZ\ C , as noted 
near the end of the previous section. For the least bound 
stars which are in the full-loss-cone regime (q(£) 1), 
Q&xi 3> 1 and a ax j w a « q. In this case, the capture rate 
does not depend on the diffusion coefficient T>{£) but 
only on the value of 1Z\ C . The boundary condition (62a) 
states that the flux T\ w is proportional to the average, 
isotropized value / ~ f\ w , and V cancels out. 

In the opposite case q(£) < 1, the capture rate is lim- 
ited by diffusion, and f(7Z) is no longer close to isotropic. 
If g < 1, a ax ; is also < 1, and the denominator in the 
expression for the capture rate (73) tends to some con- 
stant value which depends only on lZ scp , but not on q or 
1Z\ C (provided that 1Z\ C <C lZ sep ). It is largely irrelevant 
whether the boundary condition itself corresponds to the 
empty (g ax j < 1) or full-loss-wedge regime. In other 
words, in this diffusion-limited regime (both in spherical 
and axisymmetric cases) the boundary conditions (43) 
or (62a) determine f\ c for a given flux J : \ Cl which itself 
is set by the gradient of the overall steady-state profile 
of solution. By comparison, in the spherical case the de- 
nominator in the expression for the capture rate (51b) 
also depends only weakly (logarithmically) on 1Z\ C and 
is almost independent of a. Therefore, the difference 
in capture rates between the axisymmetric and spheri- 
cal problems, which results from the difference between 
ln7\L ft'(7\L sop ,7\Li c ) and ln7\Li c , is at most a factor of few 
in the case q{£) < 1. 

6. COMPARISON WITH TV-BODY SIMULATIONS 

To test the predictions of the Fokker-Planck models, 
we carried out a series of TV-body integrations of both 
spherical and flattened models of galaxies containing cen- 
tral point masses. The density model was assigned to be 
a flattened modification of the spherical Dchnen (1993) 
profile: 



p(x) = 



47T 



1 



(3- 7 ) rT(l + r) 



4—7 



l + e d 



with 7 = 3/2. This model deviates from the scale- free 
profile of equation (la) at large radii but is close to it 
inside n n fl. 

Our Fokker-Planck models are valid only for scale-free 
density profiles and at radii inside the SBH sphere of in- 
fluence. Similar iV-body studies (Brockamp et al. 2011; 
Fiestas ct al. 2012) typically assign a mass to the SBH 
particle of <~ 10~ 3 — 10~ 2 times the mass in stars, similar 
to the observed ratio. Here (as in Komossa & Mcrritt 
2008) we adopt larger values for this ratio in order to 
study in detail the region inside the influence sphere. We 
used two values for M # : 0.1 and 0.02 times the mass in 
stars (the latter set to unity). In order to simulate vari- 
ous evolutionary regimes (e.g. empty/full loss cone) we 
varied the radius r\ c at which stars are captured between 
10~ 5 and and 2 x 10~ 4 (in units of the Dehnen-model 
scale length), and we also varied the number of parti- 
cles in the system: TV = 2.5 x 10 4 , 10 5 and 2.5 x 10 5 . 
We stress that in no case would our models correspond 
to real galaxies (the capture radius is too large and the 
number of stars too small), but once we understand the 
dependence of the evolution on these parameters, we can 
scale the results to real galaxies. We summarize the pa- 
rameters of our models in Table 1. 

For the flattened models we adopted a density axis ra- 
tio of p = 0.75, which corresponds to lZ scp ~ 0.29 via 
equations (le), (lib). We did not vary this parameter 
since the foregoing analysis indicated a rather weak de- 
pendence of the flux on 7^. sep (e.g. equation 73). 

The flattened models were constructed via orbit su- 
perposition (Schwarzschild 1979) using <~ 10 5 orbits. 
While there is a unique two-integral distribution func- 
tion f(£,L z ) that self-consistently reproduces a given 
p(R,z) (Lynden-Bell 1962; Hunter & Qian 1993), there 
are infinitely many three-integral distribution functions 
(e.g. Dchnen & Gerhard 1993). In order to construct 
models that were "most similar" to isotropic spherical 
models, we chose the orbital weights in such a way as to 
minimize a global measure of the "velocity anisotropy" 
/3=1 — (<Tg + a^)/(2af). The resultant models are char- 
acterized by a non-trivial dependence of / on the third in- 
tegral, since in a two- integral, f(E, L z ) model, velocities 
are forced to be isotropic in the meridional plane only, 
i.e. <jg — of, and f3 ^ in general. Among the numer- 
ical checks that we carried out was to construct spher- 
ical models using both orbital superposition, as well as 
Eddington's inversion formula (e.g. Merritt 2013, equa- 
tion 3.47); no noticeable difference was found in the TV- 
body evolution. 

The models were evolved using the direct-summation 
TV-body code cf>GRAPEch (Harfst et al. 2008), which 
uses algorithmic regularization (Mikkola & Merritt 2006, 
2008) to increase the speed and accuracy of particle ad- 
vancement near the SBH, and includes an option for cap- 
turing particles that pass within a specified distance from 
the SBH; the mass of captured stars is added to M.. 
Integrations were carried out both on GRAPE worksta- 
tions and with the GPU-accelerated SAPPORO library 
(Gaburov et al. 2009). The accuracy parameter of the 
Hermite integrator was set to r\ — 0.01 and the grav- 
itational softening length e was set to zero. We used 
purely Newtonian gravity, as the influence of relativis- 
tic effects on the total capture rate (due mainly to stars 



15 



TABLE 1 

Parameters of the TV-body and Fokker-Planok models 

N is the number of particles, r lc is the loss cone radius (distance to SBH at which stars are captured), T re i is the central relaxation time 
defined in equation (55), T s i m is the duration of the simulation, InA is the Coulomb logarithm, r m is the influence radius, f"iocal (global) arc 
radii corresponding to energy £i OC ai( g iobal) °f transition between empty and full-loss-cone regimes for the spherical problem defined at the 
end of §4.1, JV ca pt is the total number of particles captured by the end of run (separately for spherical and axisymmetric case with axis 
ratio of 0.75, and for Fokker-Planck and iV-body models). 

Ncapt, spherical N cap t, flattened 

Model N M. ric T rel T sim InA r m nocal(global) F-P ^V-body F-P jV-body 
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Fig. 7. — Comparison of theoretical diffusion coefficients with 
results from the ./V-body simulations (spherical Dehnen model with 
M. = 0.1, N = 10 5 ). Theoretical coefficients are shown by the 
solid (7?,) and dashed (£) curves. The diffusion coefficient in energy 
has been shifted downward by two decades to avoid overlap. 

with a ~ r m ) is likely to be negligible (Appendix A). 
The integration time T s j m in A-body units was chosen 
to be a substantial fraction of the central relaxation 
time T ro i (55), but not longer, in order to avoid sig- 
nificant changes in the density profile. In all integra- 
tions, the mass of captured stars was a small fraction 
of M,, therefore we did not change r\ c as a function of 
time. A star was considered to be captured if its an- 
gular momentum near periapsis passage was less than 

used the 

latter approximation which is valid for highly eccentric 
orbits. 

Corresponding Fokker-Planck models were constructed 
using the analytically solvable equivalent Id prescription 
from the previous section, with lZ scp ,lZ\ c and a speci- 
fied at each value of £ (equations 62b, 72, 73). We used 
a number of criteria for comparing the results from the 
Fokker-Planck and A-body models. Below we present 
comparison for some of these criteria evaluated for model 
Ml (spherical and axisymmetric cases), although the re- 
sults were similar for other models. 

The first indicator is the rate of relaxation in energy 
and angular momentum. In the simulations, we sorted 



all particles in initial £ and 1Z and divided them into 100 
bins, then computed ATI 2 and A£ 2 j£ 2 for each particle 
and averaged these values within each bin. For a diffusive 
process, these quantities should grow linearly with time 
and so we fitted the time dependence with a straight line 
and took the slope of this fit as the measured value of 
(A£ 2 ,1Z 2 ). Theoretical diffusion coefficients were com- 
puted by averaging the local coefficients over the volume 
of phase space accessible at a given energy (equation E4) . 
The comparison between these theoretical coefficients 
and the values measured from the simulations is shown 
in Figure 7 for the (spherical) model Ml. The Coulomb 
logarithm is the only adjustable parameter in this com- 
parison, and the best agreement was obtained setting it 
roughly equal to the logarithm of the number of par- 
ticles inside the influence radius, InA « ln(0.3M./m*) 
(e.g. Merritt 2013, equation 5.35). Agreement was found 
to be good beyond and just inside the SBH sphere of 
influence, £ ~ 1, while at smaller radii (larger binding 
energies) relaxation in angular momentum was found to 
be faster than predicted, which may be an indication of 
resonant relaxation (Rauch & Tremaine 1996), although 
the diffusion rate was not as high as measured by Eilon, 
Kupi & Alexander (2009). For a power- law cusp with 
7 = 3/2, both ATI 2 and A£ 2 /£ 2 should tend to a con- 
stant limit for £ — > oo. Since the number of stars in the 
simulation is finite and there is a maximum value of the 
binding energy, £ « 100, for an JV = 10 5 particle sys- 
tem, we introduced an upper energy cutoff in the distri- 
bution function in computing the theoretical coefficients 
(cf. Bar-Or et al. 2012, Appendix D), which results in the 
decline of the coefficients at large £. We did not attempt 
to study the resonant relaxation in more detail since it is 
not well described by our Fokker-Planck formalism and 
more sophisticated statistical models may be needed (e.g. 
Madigan et al. 2011). In any case, the enhancement in 
the capture rate due to resonant relaxation is small due 
to the small number of particles at high binding energies 
(Hopman & Alexander 2006). 

Next we compare the properties of captured particles 
and the population of the loss cone with the predictions 
of the Fokker-Planck models, in both spherical and ax- 
isymmetric models. For every captured star, we recorded 
the energy and angular momentum at the moment of 
capture, then looked back to find their changes since the 
previous periapsis passage. Figure 8 plots changes in 
squared angular momentum during the final orbit ver- 
sus particle energy and compares it with the expected 
(average) change due to diffusion. In the case of spheri- 
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Change in squared angular momentum during the final orbit before capture, for spherical (left) and axisymmetric (right) models 



Fig. 

with M. = 0.1, N = 10 b , r lc 



10- 



Solid lines are predictions from loss cone theory: equation (74a) for the spherical case, equation (76) 



for the axisymmetric case; the latter should be regarded as an upper limit for the reasons discussed in the text, 
the loss cone boundary, L 2 = 2GM,r\ c 

cal models (Figure 8a), those changes were predicted in 
terms of q(£ ) (42) as follows: 



Horizontal dashed line is 



Aft = y/VRRTrsd = J (ftlc + Aft) T rad = 



= Vgftic(fti7+ Aft) = gfti 



l + vO + 4^ 



AL 2 = L 2 circ ATZ = 2GM.n c q 



1 + v/l + 4g- 



(74a) 
(74b) 



In spite of the substantial scatter, the measured angular 
momentum changes are well described by this approxi- 
mation. 

For the axisymmetric case (Figure 8b) the angular mo- 
mentum changes during the final orbit are higher due 
to torques from the flattened potential. We can esti- 
mate (Aft) by approximating the time evolution of ft 
(13) near the minimum as a parabola, taking ft ma x = 
and evaluating the difference 



7? 7? 

/v scpi /v mm 



Aft = (ft n 



ft • ) 



7^2 

prcc 



[(Tic + T rad ) 2 - T 2 ] , (75) 



where T\ c = (T prec / it) yj ft lc /ft max denotes the elapsed 
time after entering the loss cone until reaching the min- 
imum ft. This estimate gives an upper limit to Aft, 
expressed in terms of coefficient g a xi (59): 



Aft = fti c x 7rg axi (2g axi + n), 



(76) 



which is plotted as a solid line in Figure 8b; the measured 
values of Aft indeed lie below this upper limit but are 
higher than in the spherical case. 

The population of the loss cone in the iV-body simula- 
tions was computed as the instantaneous number of stars 
having angular momenta less than L\ c . The correspond- 
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Fig. 9. — Instantaneous number of stars in the loss cone as a 
function of energy. Solid curves are derived from the ./V-body sim- 
ulations: top (blue) - axisymmetric; bottom (magenta) - spherical 
(same model as in the previous figure). Green dot-dashed curve 
is for a full loss cone (equation 78a with AT(£) equal to its initial 
value) . Dotted red curve is the "real" , stationary loss cone popula- 
tion (equation 78b) and the dashed red curve is the time-dependent 
population, both for the spherical problem. 

ing quantity in the Fokker-Planck models is 

N <lc (£)d£= Af(£,'Jl)d'R.d£. (77) 

Jo 

We distinguish between Nf u u \ c - the number of stars if 
their distribution in squared angular momentum is uni- 
form (i.e. if the loss cone is full and the value of distribu- 
tion function Af\ c for ft < fti c is the same as the isotropic 
value AT), and 7Y roa i i c - the number of stars if Af(£,TZ) 
is taken from the true solution or its quasi-steady-state 
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Fig. 10. — Cumulative number of captured stars in the sim- 
ulations (run Ml) as a function of time (points), compared 
with predictions from the time-dependent Fokker-Planck models 
(dashed lines). Top (blue) curves: axisymmetric model; bottom 
(red/magenta) curves: spherical model. 

approximation (equation 51a). Neglecting the variation 
of orbital period with 1Z at given £, the former quantity 
becomes 

tf&U ic(£) « ndV(£ ), (78a) 
and the latter (steady-state) is 

exp(— a) — 1 + a 



iVreal lc = Mull lc 



In 



1 



(78b) 



where a(q(£)) is given by (45). Figure 9 shows the 
distribution of loss cone stars in energy in the simula- 
tions; overplotted are curves corresponding to full and 
real loss cone in time-dependent and steady-state spher- 
ical Fokker-Planck solutions. We did not derive analo- 
gous expressions for the axisymmetric case, but the re- 
sults from simulations indicate that in that latter case 
the number of stars in the loss cone is somewhat higher 
than in the spherical system for energies £ > £ global) i- e - 
above the transition from full to empty-loss-cone regimes 
for the spherical system. 

Finally, we consider the capture rate, or, rather, the cu- 
mulative number of stars captured since the beginning of 
the simulation as a function of time. The stationary so- 
lution of the Id Fokker-Planck equation underestimates 
the capture rate at early times when the phase space 
density near the loss cone has not yet reached its steady- 
state value, so we adopted the time-dependent solution 
as a basis for comparison. Figure 10 shows that the cap- 
ture rate decreases with time, i.e. the cumulative num- 
ber of captured stars grows more slowly than linearly, as 
expected. Figure 11 shows the distribution of captured 
particles in energy, which matches the Fokker-Planck so- 
lution very well, apart from an excess of captured par- 
ticles at high binding energies in the simulations, which 
is a consequence of the higher rate of diffusion in angu- 
lar momentum discussed above (Figure 7). Overall, the 
capture rate for model Ml in the axisymmetric case is 
~ 50% higher than for the spherical model with the same 



Fig. 11. — Distribution of captured stars in energy (plotted is 
the number of stars captured per unit time). Solid curves are from 
the W-body integrations: top (blue): axisymmetric; bottom (ma- 
genta): spherical models, in both cases with M, = 0.1, -/V = 10 5 , 
r\ c = 10 — 4 (model Ml). Overplotted are predictions from the 
Fokker-Planck models: dotted: stationary flux (equation (51b) for 
the spherical case (red), equation (73) for the axisymmetric case 
(blue) with the effective capture boundary 7?. e ff given by equa- 
tion (72)); dashed: time-dependent solution of the Id spherical 
problem for the same two cases. 



density profile, confirming the predictions of the Fokker- 
Planck study. The same was found to be true in the other 
models of Table 1: flattening was never found to make 
more than a factor of two difference. Moreover, models 
for which the transition to the full-loss-cone regime in 
the spherical problem occurs well within the radius of 
influence (M2, M3, M5) showed less difference than the 
model M4 which is mainly in the empty-loss-cone regime, 
which is in line with theoretical predictions at the end of 
previous section. 

We also recorded the orbital parameters of captured 
stars at their final apoapsis passage and followed the or- 
bits in the smooth potential used to construct the flat- 
tened models. Figure 12, top panel, shows that most of 
these stars found their way into the SBH while being on 
saucer orbits, a result also predicted by the axisymmetric 
Fokker-Planck models. 

As a final remark, we tested that the iV-body models 
were in dynamical equilibrium by examining the evolu- 
tion of Lagrangian radii of shells containing given frac- 
tions of the total mass (5%, 10%, etc.), and also the 
axis ratios of the models. In integrations with captures 
disabled, these did not change apart from small fluctu- 
ations. When captures were enabled, Lagrangian radii 
expanded slightly with respect to time (corresponding to 
the energy input from the SBH), while the axis ratios 
did not change appreciably. The SBH particle did not 
remain precisely at the model center but rather experi- 
enced Brownian motion (Merritt, Berczik & Laun 2007); 
however the amplitude was at least an order of magni- 
tude smaller than the influence radius. Brockamp et al. 
(2011) found no substantial differences in capture rate 
between simulations with fixed and wandering SBHs. 
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Fig. 12. — Top panel: Distribution of captured stars (per unit 
of time) by orbital type for axisymmetric model Ml. Dotted red: 
saucers; dashed green: tubes; dot-dashed blue: chaotic; solid line: 
total. 

Bottom panel: Fraction of different types of orbits in the axisym- 
metric Schwarzschild model Ml, as a function of energy. Dotted 
red: saucers; dot-dashed blue: chaotic (numbers are on the right- 
hand ordinate); solid green: sum of these two; dashed black: cu- 
mulative mass (fraction of stars with energy > S, numbers on the 
left ordinate). 



7. THE ROLE OF CHAOTIC ORBITS 

Apart from some special, fully integrable cases (e.g. 
Sridhar & Touma 1997), most axisymmetric potentials 
containing central point masses are characterized by 
chaotic motion in the low-angular-momentum parts of 
phase space beyond the influence radius (Figure 12, bot- 
tom panel). Chaotic orbits still respect two integrals of 
the motion, £ and 7Z Z , but in the absence of a third inte- 
gral, they can in principle fill the accessible region in the 
meridional plane, allowing them to be captured as long 
as 1Z Z < 1Z\ C - We estimate the capture rate from these 
orbits by the following argument, similar to an argument 
of Magorrian & Tremaine (1999). 

Assume that the chaotic orbits occupy a region in the 
1Z — 1Z Z plane with 1Z < lZ c h- The value of lZ c h plays 
the same role as !Z S ep inside the radius of influence, and 
is comparable to it for the same degree of flattening. 
Furthermore, we assume that every chaotic orbit with 
a given 1Z Z can attain values of 71 € [lZ z ..7Z sep ] with 
equal probability. Recall that the number of stars with 
given {1Z,1Z Z }, which may be identified with the proba- 
bility of finding a star in a given interval of dlZdlZ z , is 



dj\f = Gsph f(7Z, TZ-z) dlZ dlZ Zl with G sp h given by equa- 
tion (31). Then the fraction of time such an orbit spends 

below the capture boundary is (s/1Z\ c — ^J1Z Z ) / ^f!Z c \i, and 
this is essentially the probability of being captured dur- 
ing one radial period (assuming the full-loss-cone bound- 
ary condition, i.e. that the change in angular momen- 
tum during one period is much larger than the capture 
boundary, which is reasonable for chaotic orbits). 

Next we evaluate the time-dependent rate of capture 
of stars from chaotic orbits. From the above argument 
it follows that we need to consider f(lZ z ) decaying expo- 
nentially at every value of 1Z Z from its initial value /i n it , 
but with a different rate: 



f{U z ,t) = /i„it exp 



t 



rati 



The total number of chaotic orbits with 1Z Z < 1Z\ C and 
their capture rate is then given by 

N ch (s , t)d£= J ' c g £ f(n z ,t) (yn~Jn z - 1) dn z ds 
* Qsfimt iVwiTc 1 z cx p(- 2t ) dS> (79a) 

F*{8, t) dE = G£f »^ 1-^ + 1)cM-2t) d£ 



rad 



2r 2 



(79b) 



where 



7~ — ^/^drain ) ^drain = 2T ra d 

V^ch/^lc- (80) 

From here it is clear that if we identify /i lut with the ini- 
tial value of the distribution function in the loss cone f\ c , 
then the capture rate is initially equal to the draining rate 
of a uniformly populated loss cone (48). In particular, 
when /i n it = /, we recover the standard, full-loss-cone 
draining rate, regardless of the value of 7?. c h. On the 
other hand, the draining time does depend on 7?. c h, since 
the number of stars in the chaotic region to be drained is 
2-\/7?. c h/7?-ic times larger than the number of stars in the 
loss cone, therefore the draining time is longer than the 
radial period by the same factor. At times much longer 
than the draining time, the capture rate declines as t~ 2 , 
not exponentially, since it is dominated by the draining 
of chaotic orbits with 1Z\ C — 1Z Z 1Z\ C . 

There is a great deal of similarity between the cap- 
ture rates from the loss wedge of the saucer region of 
phase space for regular orbits within the radius of influ- 
ence, and the chaotic region outside it. The details of 
draining are somewhat different (in particular, for the 
regular orbits the draining rate declines as i~ 3 , as noted 
by Magorrian & Tremaine (1999)), but since the draining 
time for saucer orbits is usually much shorter than Hub- 
ble time, we ignore that distinction and adopt the same 
expressions as for chaotic orbits. In both cases, the lo- 
cal boundary condition for the loss region corresponds to 
the full loss cone (62a), at least for the case g ax i > 1 rele- 
vant for all but the most tightly bound orbits. However, 
the global shape of the steady-state solution depends on 
whether the overall flux into the low angular momentum 
region is limited by diffusion {q{£ ) < 1) or not. In the 
first case, the steady-state solution will still have a log- 
arithmic form for 1Z > lZ c h or lZ SC p, corresponding to 
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TABLE 2 

Comparison of three geometries 

Spherical Axisymmetric Triaxial 

Fraction of stars 

with 7?. m i n < 1Z\ C 1Z\ C V^lc^-sep Tlsep 

Draining time T drain T rad > T prcc > Tprcc 

some effective capture boundary TZ e s, and the capture 
rate depends on this effective boundary only logarithmi- 
cally. In the latter case, the capture rate is essentially 
the full loss cone rate for an isotropic distribution func- 
tion. The latter case, however, is rarely attained because 
q(£) rapidly drops with decreasing binding energy On 
the other hand, if the draining time for chaotic orbits is 
comparable to the Hubble time, then their capture rate 
may still be quite high even in the absence of relaxation, 
provided that the initial value /i n it of the distribution 
function inside the chaotic region was not much different 
from the isotropic value /. 

The combined effect of draining and relaxation may 
be approximately accounted for by the following recipe. 
Let J" re i(£, t) d£ be the capture rate per unit energy from 
the Fokker-Planck equation with initial conditions corre- 
sponding to the loss region being initially empty (i.e. the 
solution considered in § 5). Since the loss region initially 
may have some nonzero value of /, < /i n ; t < /, the 
phase-space gradient of / near the loss region bound- 
ary will be less than arising in the Fokker-Planck so- 
lution, and the capture rate from relaxation alone may 
also be lower. We approximate the total capture rate by 
the sum of the draining rate J-drain and the collisional 
flux ,F re ] multiplied by 1 — -/V c h(i)/-/V c h.o, where is 
the number of chaotic orbits at t = with /j n j t = /. 
This expression is used to compare Fokker-Planck mod- 
els against ./V-body simulations in §6 and to compute the 
capture rates for real galaxies in §9. It is important to 
note that the effective capture boundary TZ c s defined in 
(72) , as the parameter controlling the overall shape of the 
steady-state solution and the gradient of the distribution 
function (and hence the capture rate due to diffusion), 
is not the same as the size of the loss region ( ^/7^ sop 7^i c 
for saucer orbits, 2yJTZ c \JZ\ c for chaotic orbits) which de- 
termines the draining time. The former, being a fixed 
fraction of 7?. sep (R- C h), is usually much larger than the 
latter. 

8. TRIAXIALITY 

For completeness, and to put our results in a broader 
context, we briefly discuss the case when the stellar cusp 
around the SBH is triaxial. Triaxial potentials support 
two distinct familes of tube orbits, circulating about the 
long and short axes of the triaxial figure. In addition, 
there is a new class of centrophilic regular orbits, the 
pyramids (Merritt & Valluri 1999, Figure 11). The defin- 
ing feature of pyramids is that lZ m i n = for all of them, 5 
and a star on such an orbit will eventually get into the 
SBH even without the assistance of collisional relaxation. 
The fraction of phase space occupied by pyramids is com- 
parable to that of saucer orbits, i.e. ~ 7^. sop . Outside 
the radius of influence, the regular pyramid orbits are 

5 Rclativistic precession alters this conclusion for the most bound 
orbits (Merritt & Vasilicv 2011). 



replaced by chaotic orbits, which are however still cen- 
trophilic (Poon & Merritt 2001). 

Table 2 summarizes the three cases. The number of 
stars that potentially can be captured (loss region, stars 
with 7?. m i n < 1Z\ C ) increases with decreasing symmetry, 
however the instantaneous number of stars in the loss 
cone is the same (the fraction of time a star on a loss 
region orbit actually has instantaneous 1Z < 1Z\ C exactly 
balances that). Consequently, the survival time of these 
stars also increases with decreasing symmetry: Td ra i n is 
larger than T ra( j by the same factor as the number of 
stars in the loss region versus the loss cone, assuming a 
full-loss-cone draining rate. 

In the absence of relaxation, the loss region is rapidly 
depleted in both spherical and axisymmetric cases (the 
latter - except for the most massive SBHs), but in the 
triaxial case the draining time of the loss region may 
be comparable to or even exceed galaxy lifetimes (Mer- 
ritt & Vasiliev (2011); in that paper a rather small de- 
parture from spherical symmetry was considered; for 
~ 0.1 or for more massive SBHs the time will be 
longer.) If, as suggested by Merritt & Poon (2004) and 
Holley-Bockelmann & Sigurdsson (2006), the capture of 
stars from centrophilic orbits (both regular pyramids in- 
side the radius of influence, and chaotic orbits outside 
it) can sustain a full-loss-cone feeding rate for an ini- 
tially isotropic distribution of stars, then it will remain 
near this level for a time ~ Td ra i n - For stellar systems 
older than that, it is necessary to take 2-body relaxation 
into account, and the outcome will probably be similar 
to what happens in the g a xi > 1 , Q < 1 regime of the ax- 
isymmetric problem: reshuffling of stars in angular mo- 
mentum near the loss cone boundary due to nonsphcrical 
torques is efficient enough to keep the loss cone full, but 
the value of f\ c near loss cone is much smaller than the 
average (isotropized) value /, because the supply of stars 
into the \cw-lZ region is limited by the diffusion from 
higher 1Z. Extrapolating the estimates of draining times 
for axisymmetric galactic models from the next section 
into the triaxial case, one may conclude that for most 
massive galaxies the lifetime of centrophilic orbits may 
indeed be longer than Hubble time, provided that the 
triaxiality is not destroyed by effects of chaos (Merritt & 
Quinlan 1998). 

We stress that genuinely centrophilic orbits exist only 
in the triaxial case, since in the axisymmetric geometry 
the conservation of L z precludes orbits from reaching 
arbitrarily small radii. However, some degree of non- 
axisymmetry is to be expected in every real galaxy. 

9. ESTIMATES FOR REAL GALAXIES 

Having verified that the Fokker-Planck models agree 
well with the iV-body simulations in the range of pa- 
rameters feasible for the latter, we now extrapolate the 
Fokker-Planck theory to parameters characteristic of real 
galaxies. To that end, we consider a family of models rep- 
resented by Dehnen (1993) density profiles having inner 
cusp slopes 0.5 < 7 < 2, and take the SBH mass to be 
10~ 3 of the total galaxy mass. The corresponding in- 
fluence radii lie well inside the break radius separating 
the inner p ~ r~ 7 cusp from the outer p ~ r~ 4 profile, 
so that only the density normalization at the radius of 
influence (say) matters. That density is, in principle, 
an independent parameter but we fix it via the require- 
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Fig. 13. — Estimates of capture rates M depending on the black hole mass Mm and the parameters of galaxy. 
Left panel: Predictions from Fokker-Planck models for galaxies with inner density cusp slope 7 = 1.5 for M, < 10 s Mq and 7 = 1 for 
M, > 10 8 Mq (the discontinuity corresponds to ~ 3 times higher rates for more concentrated galaxies). Galaxy models are scaled to match 
M. — <t relation (81) with a = 8, = 4.5. Solid red curve and dashed purple lines: stationary and time-dependent flux in the spherical 
case; dot-dashed and dotted blue lines: same for axisymmetric case with TZ Bcp = 0.1; double-dashed blue line: flux in the axisymmetric 
case with contribution from loss region draining. On average, stationary capture rates differ by a factor of 2-3 between spherical and 
axisymmetric cases; time-dependent solution generally gives higher rates since it starts from an initial condition with strong gradients in 
the distribution function near the loss region, which may not be physically motivated and hence provides an upper limit to possible rates 
in realistic conditions. 

Right panel: Comparison of stationary capture rates in the spherical case for varying density cusp slope, with shaded regions corresponding 
to uncertainties associated with different parameters a, f} in M, — a relation (81), Table 3). From top to bottom -7 = 2, 1.5, 1, 0.5. Points 
are taken from previous studies: purple open boxes — Syer & Ulmer (1999); blue upward and downward triangles — Magorrian & Tremaine 
(1999) for the spherical and axisymmetric cases; black open and filled circles — Wang & Merritt (2004) for cored and cuspy galaxies; black 
double-dashed line - their analytic estimate for a singular isothermal sphere (7 = 2); red dot-dashed line - Brockamp et al. (2011) using 
their extrapolation from TV-body simulations of n = 4 Sersic galaxy models; green crosses - Fiestas et al. (2012) for TV-body evolution of 
King models scaled to Milky Way nuclear star cluster with a 7 = 1.75 cusp (spherical and flattened case). 



TABLE 3 
Parameters of M, — a relation 



a f3 Reference 

8.11 4.72 Merritt & Ferrarese (2001) 
8.13 4.02 Tremaine et al. (2002) 

8.12 4.24 Giiltckin et al. (2009) 

8.13 5.13 Graham et al. (2011) 
8.29 5.12 McConncll et al. (2011) 



ment that the galaxy satisfy the so-called M. — a relation 
linking M m to the velocity dispersion observed near the 
center of the galaxy. The M, — a relation is typically 
written in the form 

log A/. = a + /31og(cr/200km s" 1 ), (81) 

with parameters a~8,4</3<5, depending on the 
study. An incomplete selection of papers during the last 
decade is summarized below: 

We scale the velocity unit of the Dehnen model by 
identifying the line-of-sight velocity dispersion, a p (R), at 
the influence radius with the quantity a in equation (81); 
this is reasonable given that a p is a weak function of pro- 
jected radius R for r m < R < R e , with R e the half-light 
radius. (In the case 7 = 2 the central a p in the model 
lacking a SBH was used.) Our models are thus defined 
by the two parameters {M,,j}. Setting a = 8,/3 = 4.5, 
we have r m = {50,35,18,6} x (M./10 8 M Q ) 56 pc for 
7 = {1/2, 1,3/2, 2}. 

The radius r\ c that defines the loss sphere around the 



SBH is the larger of the radius of tidal disruption, r t id, 
and the (Newtonian) periapsis of an orbit that just con- 
tinues inside the event horizon. For the eccentric orbits 
that dominate the flux into the SBH, the latter quantity 
is ~ 8r g = 8GM,/c 2 for a Schwarzschild (nonrotating) 
SBH (Will 2012); this is the radius of periapsis of a Ke- 
plerian orbit having the critical angular momentum. A 
star is disrupted if the periapsis radius is less than 

nn 2/3 ( M. V 2/3 fm*\~ 1/3 r* 

(82) 

Here r\ depends on the stellar equation of state and is 
~ 0.84 for a solar-type main-sequence star. These dis- 
ruption events, as opposed to direct captures, may be 
observed as optical and x-ray flares in otherwise quies- 
cent galactic nuclei (Strubbe & Quataert 2009). From 
the condition r t id > 8r g we find that solar-type stars on 
eccentric orbits are disrupted (not swallowed) if M. < 
1.2 x 10 7 M©; disruption can occur for any M. < 10 8 A/ Q 
if the star is on a less eccentric orbit, or for Kerr SBHs 
even more massive than 10 8 M Q (Kesden 2012). Red 
giants or AGB stars can also be disrupted (or at least 
tidally limited) by SBHs more massive than 10 8 Mq 
(MacLeod et al. 2012). In what follows, we compute 
the total number of events associated with a given n c ; 
the ratio of number of tidal disruption flares to the total 
number of capture events is well studied in the literature 
and we do not consider it separately here. 
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We used Fokker-Planck models to evaluate steady- 
state and time-dependent capture rates M for galaxies 
after 10 10 years starting from an initially isotropic distri- 
bution function, in both the spherical and axisymmetric 
geometries (using our one-dimensional approximation of 
§5), for lZ scp =0.1. The latter value is meant to represent 
a "typical," moderately-flattened system; the results do 
not strongly depend on 1Z SC p- In the time-dependent cal- 
culations, the initial conditions consisted of the isotropic 
models with loss-cone orbits removed; as a result, these 
initial models are characterized by strong gradients of / 
with respect to L near the loss cone. 

The left panel of Figure 13 shows results for two series 
of models: models with a steep (7 = 1.5) central cusp and 
M, < 1O 8 M ; and models with shallow (7 = 1) cores 
and M, > 1O 8 M . Over the entire range of M # , the 
steady-state capture rates differ by only a factor of 2 — 3 
between spherical and axisymmetric geometries, consis- 
tent with the discussion near the end of §5. This result 
holds for any galactic model and depends only weakly 
on lZ sep . The time-dependent rates are generally higher 
than in the steady state, due to the strong gradients 
in the initial conditions. Especially for massive galax- 
ies with long relaxation times, the approach to a steady 
state is slow and the flux at early stages is much higher 
than in equilibrium. 

For the most massive SBHs (M. > 10 9 Mq) the drain- 
ing time of the loss region becomes comparable to the 
Hubble time and the capture rate is dominated by drain- 
ing of chaotic orbits (double-dashed line) , reaching values 
up to 10~ 3 M Q yr _1 at the upper end of the M. range. 
However, this should be regarded as a strong upper limit 
since we do not know the initial state: for instance, if the 
SBH formed as the result of a merger of a binary SBH, it 
is very likely that the low angular momentum region of 
phase space will have been depleted in the course of the 
binary's evolution, and the capture rates could initially 
be much lower than the steady-state values (Milosavl- 
jevic & Merritt 2003; Merritt & Wang 2005). 

The dependence of M on M. for massive SBHs can 
be simply estimated as follows (e.g. Syer & Ulmer 1999). 
The nuclei of massive galaxies are in the empty-loss-cone 
regime, so the flux per unit energy is roughly F{£) d£ ~ 
N{£) d£j [T re i(£) \n"R^}{£ )] . The capture rate peaks at 
r w r m (Merritt 2013, section 6.1.4.1) so the total flux is 
estimated as M ~ ^ r (finfl)^infi ~ M,/(T rc \ InTvLjT, 1 ) ~ 
(m*/M.) a 3 /G. Plugging in the M, — a relation we 

obtain N oc M. 3//3_1 , and even the normalization con- 
stant evaluates to a reasonable ~ 10~ 5 M Q yr _1 for 
M, = 10 8 Mq, despite the crudeness of the estimate. 

The right panel of Figure 13 shows the uncertainties in 
the capture rate associated with the parameters a and 
(3 in the M. — a relation (81) and the slope of the den- 
sity cusp 7. Plotted are stationary capture rates for the 
spherical case; other values scale roughly in the same 
proportion. The capture rates evaluated for a selection 
of individual galaxies from several previous studies are 
also plotted for comparison. It is clear that the scatter 
in the derived values is fairly large, about two orders of 
magnitude, although a general trend of decreasing rate 
with increasing M, is clear. Comparing the inverted tri- 
angles in Figure 13b with the double-dashed curve in 



Figure 13a, we see that our estimates for the capture 
rate due to draining of chaotic orbits are substantially 
higher than those of Magorrian & Tremaine (1999). As 
argued above, this is most likely an overestimate result- 
ing from simplistic initial conditions. It is also due partly 
due to our selection of a different relation between the 
capture rate and the isotropized distribution function; 
had we used their expression (66) instead of our (71), 
the steady-state flux in the axisymmetric case would be 
factor of a few lower, although it should not affect the 
draining rate of centrophilic orbits which starts to dom- 
inate the capture rate at M. > 10 9 Mq. 

Overall, it is fair to say that our estimates predict cap- 
ture rates in the range 10 -5 — 10~ 4 Mq yr _1 for less mas- 
sive galaxies, and a fewx 10~ 6 — 10~ 5 Mq yr _1 for giant 
galaxies with SBH masses in excess of 10 8 M Q . This is 
roughly consistent with the observationally derived es- 
timates of rates of tidal disruption flares (Donley et al. 
2002; Gezari et al. 2008; van Velzen & Farrar 2012). (Ax- 
isymmetric) nuclear flattening may increase these num- 
bers by a factor of few, and triaxiality may have a more 
dramatic effect on the consumption rate of the most mas- 
sive SBHs provided there are enough stars on centrophilic 
orbits. 

One potentially important feature of axisymmetric 
(and triaxial) systems is that most stars are consumed 
in the full-loss-cone regime of boundary conditions (in 
spite of the fact that the angular momentum need not 
be isotropic). This means that stars approach the SBH 
with a wide distribution in periapsis radii, as opposed 
to "barely touching" the disruption sphere in the empty- 
loss-cone regime. As a consequence, many stars will be 
strongly tidally distorted before disruption, which may 
result in a distinct observational signature (Strubbe & 
Quataert 2009), although more recent studies show that 
the difference may not be so pronounced (Stone et al. 
2012). 

10. CONCLUSIONS 

We have considered collisional (gravitational two- 
body) relaxation processes near supermassive black holes 
(SBHs) in spherical and axisymmetric models of galactic 
nuclei. We derived a Fokker-Planck formalism and com- 
pared its predictions with direct A-body simulations of 
capture. 

Inside the SBH radius of influence, the unperturbed 
motion of stars is regular, admitting three integrals of 
motion (energy, ^-component of angular momentum L z 
and the secular Hamiltonian %). There are two fami- 
lies of orbits, tubes and saucers; the latter exhibit large 
angular momentum variations and stars on saucer or- 
bits approach much more closely to the SBH than would 
be expected based on their average angular momentum. 
Regularity of the motion allowed us to write down the 
orbit-averaged Fokker-Planck equation and to calculate 
the diffusion coefficients based on the standard formal- 
ism. We discussed the appropriate boundary conditions 
for capture by the SBH, and numerically solved the 
two-dimensional (H — L z ) Fokker-Planck equation. We 
showed that its solution can be well approximated by an 
equivalent one-dimensional solution for diffusion in angu- 
lar momentum, with some effective boundary conditions. 
An important difference with the spherical case is that 
the boundary condition at the loss region typically cor- 
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responds to the full-loss-cone regime, in the sense that 
a change in angular momentum per one orbital period 
is larger than the size of the loss cone. Nevertheless, 
the global shape of the solution, and, consequently, the 
capture rate of stars, is determined by the diffusion co- 
efficient (inversely proportional to the relaxation time), 
and depends only weakly on the effective size of the loss 
region. This treatment was not entirely self-consistent, 
as it did not account for resonant relaxation or for the 
breakdown of the orbit-averaged approximation near or 
beyond the radius of influence, but it nevertheless sug- 
gests a conclusion which turns out to be robust: com- 
pared with the spherical case having the same real (not 
effective) capture boundary, the flux into the loss cone 
is higher in the axisymmetric case, but not by a large 
factor, and only in the regime where the loss cone would 
be empty in the spherical system. In the axisymmetric 
case, most of the stars find their way into the SBH while 
on saucer orbits. 

We also carried out a number of iV-body integrations 
which were found to agree remarkably well with the cor- 
responding Fokker-Planck models. The agreement was 
not limited to the number of captured stars; other quan- 
tities like the distribution of energy of the captured stars, 
the loss cone population, the diffusion coefficients, and 
the change in angular momentum during one radial pe- 
riod prior to capture were also found to be well repro- 
duced. We note, however, that the angular momentum 
diffusion is larger in TV-body models for stars with high 
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binding energies (close to the SBH), which may be an in- 
dication of resonant relaxation not accounted for in the 
Fokker-Planck models; nevertheless, its influence on the 
capture rate is likely to be small. 

We applied our Fokker-Planck formalism to realistic 
galaxies with SBH masses M. = (10 6 - 1O 1O )M . We 
found that stationary capture rates are in the range 
(10~ 4 — 1O~ 6 )M0 yr _1 in spherical galaxies and a fac- 
tor 2 — 3 higher in flattened systems, with an overall 
trend of decreasing event rate with increasing M.. Time- 
dependent solutions were found to give generally higher 
estimated capture rates; however, that result is likely to 
depend strongly on the assumed initial conditions. In 
particular, for massive (M, > 10 9 Mq) SBHs in axisym- 
metric nuclei, the draining time of chaotic orbits just 
outside the radius of influence can be comparable with 
the Hubble time; if such orbits were not depleted by 
(for instance) the binary SBH that preceded the single 
SBH, capture rates might reach 10 -3 Af©yr . Such high 
rates may be more relevant to triaxial galaxies in which 
the number of ccntrophilic orbits can be large enough to 
maintain a full-loss-cone capture rate for a Hubble time. 
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APPENDIX 

A. EFFECT OF RELATIVITY ON THE UNPERTURBED MODTION 

We discuss briefly the character of the motion in axisymmetric nuclei when the lowest-order post-Newtonian (PN) 
corrections are included in the equations of motion. The 1PN accelerations imply an orbit-averaged rate of of periapsis 
advance 



cku \ 
~dt ) 



GR 



2tt 3GM. 

T rad c 2 aR 



(Al) 



(Merritt 2013, equation (4.205)). The other elements of the osculating orbit exibit no secular variations at this PN 
order, and we ignore rotation of the SBH. Expressed in terms of the dimensionless time variable r = 2-Kt/Tyi defined 
in equation (9), the relativistic precession rate becomes 
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which can be added to the orbit-averaged equation of motion for lu, equation (9). Setting 7 = 1 (Q = 3/2) in the 
Hamiltonian (8), the two nontrivial equations of motion become 
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If there is a fixed point, its angular momentum can be found by setting lu = when lu = tt/2, or 
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As k is increased from zero, the value of I at the fixed point increases, reaching t = 1 when 
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For k > K2 all orbits are tubes. For k\ < k < k 2 , where 
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there is one family of tubes, passing below the fixed point in the (72., ui) plane, and a family of saucer-like orbits that 
librate around the fixed point. For k < m, the separatrix encloses the fixed point and there are two families of tubes, 
at low and high angular momenta. In the latter case, the angular momentum associated with the separatrix at lu = 
can be found by setting lu = in equation (A3b): 
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Since n\ < 1 < k 2 , saucer- like orbits are present only when 
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Even when saucers are present, their angular momentum variations are limited by the relativistic term. A rough 
lower limit on the attainable angular momentum for saucers can be derived by equating the change in i due to the 
torques over one GR precessional period to I: 
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Fig. 14. — Solutions to the equations of motion (A3) for TZ Z = 0.01, k = 0.02 and e p = 0.035. The fixed point, equation (A4), is shown 
by the cross. There are two familes of tube orbits: tube orbits above the fixed point, which are present even when n = (no GR); and 
tube orbits below the fixed point, having sufficiently small i that relativistic precession quenches the effects of torques due to the flattened 
potential. Equation (A10) predicts TL ~ 1.4 X 10~ 2 for the minimum TL reached by saucers; the actual value is ~ 2.4 X 10~ 2 . 



which yields 



c min ~ 

67T € v 



(A10) 



Figure 14 shows numerical solutions of the equations of motion (A3) for 1Z Z = 0.01 and k = 0.02; the nuclear 
flattening parameter has the same value as in Figure 2. The neglect of GR precession on the evolution of saucer orbits 
is justified if 



2n c 
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where the latter inequality expresses the fact that r\ c = 8r g for direct captures and larger than that for tidal disruptions. 
Requiring that £ min < 4y r g /a and using equations (3, A2, A10), we obtain 
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where the second line approximates r m w GM,/a 2 . This is essentially the same condition that was obtained in Merritt 
& Vasiliev (2011) for capture of pyramid orbits. Since the quantity in the brackets is likely to be small (unless e p 
is tiny), and because most of the flux into the SBH comes from orbits with a rss r m , equation (A12) suggests that 
relativity is not likely to be important for the total capture rate. 

B. ORBIT-AVERAGED HAMILTON! AN 

Here we give the exact expression for the orbit-averaged Hamiltonian corresponding to the potential of equation (1). 
It is convenient to express r and z in terms of eccentric anomaly r\ rather than mean anomaly w: 



r = a(l - e cos 77) , z = a\jl — ^| sin w(cos 77 - e) + \J \ - e 2 cos u> sin 77 , e = \J \ - P. 
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C. RELATION TO THE LIDOV-KOZAI PROBLEM 

Motion in the hierarchical three-body problem is often derived from a doubly-averaged Hamiltonian after expressing 
the equations of motion in Jacobi coordinates and retaining only the lowest-order (quadrupole) term in the perturbation 
potential. The "inner restricted problem" (Kozai 1962; Lidov 1962) assumes furthermore that the test mass orbits 
well inside the perturber mass. The averaged Hamiltonian describing the test particle, with an appropriately chosen 
unit of time, is 



5 72 ( 72 z 
H * = - 6+2 + {'-it 



5 / t~* \ 2 72 
-(l-72)sm 2 W + - 



(Cl) 



(Merritt 2013, equation 4.315) which may be obtained from equation 8b by setting Q — 5/2 and eliminating unity in 
the first bracket of the first term (equivalent to taking the limit e p — > oo, eliminating terms that do not contain e p 
and normalizing the unit of time to e p ). Equation (lib) then yields a value of 5/3 for the parameter 72 sop . As in the 
axisymmetric problem, motion in the Lidov-Kozai problem also exhibits two regimes: circulation in lo (corresponding 
to tube orbits) or libration about lo = tt/2 (analogous to saucers), the latter appearing for 72 z < 3/5. The equation of 
motion for 72 (12a) is 

(972 

-j— = -2V6 y/(K! -K)(K-K 2 )(n 3 -K) , (C2) 
ar 

with the same relation between 72-1,72.2,72-3 as in equation (12b): 



72i, 2 = 72*± V72 2 -(5/3)72 z (C3a) 
72 3 = 7fK + 5/6 + 72 z /2 (C3b) 
72* = (5 + 572 z - 272 3 )/6. (C3c) 

The case of circulation corresponds to 72-2 < 72 < 723 < 72. i, and libration to 722 < 72 < 72 1 < 723; in the latter case the 
relation between minimum and maximum values of 72. is the same as in equation (16), from which it follows that the 
librating regime exists for 72 z < l/72 scp = 3/5. The separatrix between "tubes" and "saucers" is at Hk = 1/6 — TZ z /2 
and the fixed-point saucer has Hk = 5/3 + 272 z — y/151Z z . For comparison, in the oblate axisymmetric potential 
considered throughout this paper, the first regime (tube orbits) has 722 < 723 < 72 < 72-i and the second (saucers) has 
72 3 < 72 2 < 72 < 72.1. 

D. SAUCER ORBITS BEYOND THE SPHERE OF INFLUENCE 

Beyond the influence sphere, orbits similar to the saucers can still exist in axisymmetric potentials (Richstone 1982; 
Lees & Schwarzschild 1992; Evans 1994), but they are not describable in terms of osculating Keplerian elements. 
Typically such orbits are described as tube orbits that lie close to a resonance between the radial and vertical motions. 
To make the correspondence with our work more clear, we recast the fixed-point saucer orbit near a SBH in terms of 
the Cartesian variables (R,z), i.e. cylindrical coordinates in the meridional plane. Setting lo = tt/2 in equation (Al) 
yields for the fixed-point orbit that generates the saucers 

*_ (1 .^"_ML^l».( 1 .4 T )4 ( D1 ) 

a z sin i a \ sin i J a z 

where it is understood that e and cosi have their fixed-point values: 



e 2 = 1 - y'KscpKz, cos 2 i = 



I n z 

72 S cp 



Equation (Dl) is a hyperbola in the (R, z) plane, between the points 

7?max = (1 + e)acosi, z max = (1 + e)asini, 
7?min = (1 — e)acosi, z max = (1 — e)asini. 



The curve crosses the equatorial plane at R = a^/72 S op72 z = a(l — e 2 ), the semi-latus rectum. 

When the fixed-point saucer orbit first appears, at 72 z = 72 sop , it lies in the equatorial plane (cosi = l,z max = 0) 
and has zero vertical thickness. For 72. z values smaller than this maximum, the trajectory (Dl) can be interpreted as 
a 1 : 1 resonance between motions in the R and z directions. 

Described in this way, saucer orbits near a SBH are seen to have very similar properties to orbits described by other 
authors in more general axisymmetric potentials. Lees & Schwarzschild (1992) studied orbits in scale- free axisymmetric 
models with logarithmic potentials, p ~ r~ 2 and no central SBH. For models with density axis ratio 0.265, they found 
that saucers first appear at 72 z w 0.48; for smaller 72. z the fixed-point orbit (which they called a "reflected banana") 
traces a path in the (i?, z) plane similar to equation (Dl). They noted that motion near the fixed-point orbit is regular, 
i.e. non-chaotic. Similar orbits were described by Richstone (1982) (who called them "pipe orbits") and Evans (1994) 
in surveys of orbits in other scale-free families of oblate models. 



26 



E. LOCAL DIFFUSION COEFFICIENTS 

Here we present the local (position-dependent) diffusion coefficients appearing in equation (22), expressed in the 
spherical coordinates (r,9,<fi) and generalized velocities (S = <f>(r) — v 2 /2, 1Z = kL 2 , 1Z z = kL 2 ). Wc denote 
k = L- 2 C (E), k' = dn/dE = -dn/dS, Q = 1 + (v 2 /2)(k' / 7s). 

<A£> = -«<A«||> - i(A«f) - I<A«i), (El) 

?2\ _„,2/ a„,2\ 
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(AKAK Z ) = AKK Z Q Z — ^ + 2"K"K 

ir 

These coefficients have been expressed in terms of velocity diffusion coefficients for {v\\, v±} via integrals of the distri- 
bution function /(£/) describing the field stars (of mass m*) using the relative potential * = — $: 

v{Av\\) = -(l + ^)l 1/2 , (E2) 

(A V f) = |(/o + J 3/ 2), 

(A«i) = |(2/ + 3/ 1/2 -/ 3/2 ), 

where 

7 = 167T 2 G 2 ml In A / d£ f f{£ f ), (E3) 
Jo 

J n/2 = 16n 2 G 2 ml In A jf ^— ^ J . 

When making comparisons with the relaxation rates measured from Af-body simulations we averaged the diffusion 
coefficients over the subspace £ — const: 

/•r max (f) 
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F. COORDINATES IN 2D FOKKER-PLANCK EQUATION 

Although we discuss the two-dimensional Fokker-Planck equation in the % — 1Z Z plane throughout the paper, it is 
more convenient to obtain the numerical solution using a different set of coordinates, for which we have implemented 
two variants: {y, ji) and (u, £), where 



v = H + ll z ; M= - U-J&2_ l b=— + 1 - W ; g = U = — . (Fl) 

^ V V '<-sep / 'Vscp rt — rto 1 — /v sop 
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The lines of constant v are diagonal lines in the % — 1Z Z plane, parallel to the capture boundary in the tube (Ji > 0) 
region, and lines of constant fi or £ are also straight lines in that plane, designed in such a way that the capture 
boundary in the saucer region (H > 0) has fj, = const = lZi c /lZ sep or £ = const = —TZic/Ho- This facilitates setting 
boundary conditions on these capture boundaries which are parallel to coordinate axes. We used both coordinate sets 
to cross-check the results. 



